Differential Algebraic Equations
1. Introduction to DAEs
Differential-Algebraic Equations are a class of functional equations that contain both differential equations (describing the evolution of the system) and algebraic constraints (restricting the state space). Unlike standard Ordinary Differential Equations (ODEs), DAEs are not explicitly solved for all derivatives.
A general DAE system is expressed in the implicit form: \(F(t, y, y') = 0\)
If the Jacobian \(\frac{\partial F}{\partial y'}\) is non-singular, the system is essentially an implicit ODE. If it is singular, the system is a “true” DAE.
2. The Concept of Index
The difficulty of solving a DAE is measured by its index. The most common definition is the differentiation index: the number of times you must differentiate the algebraic constraints to express the system as a set of explicit ODEs. * Index 0: An ODE. * Index 1: The most common solvable DAE (e.g., the algebraic variables can be solved for directly). * Higher Index (2+): These are numerically unstable and usually require index reduction techniques before solving.
3. Solving DAEs with sepalsolver
In modern computational environments like C#, DAEs can be solved using the SepalSolver library, which utilizes a Mass Matrix formulation: \(M y' = f(t, y)\) Where \(M\) is a singular matrix.
—
4. Examples and Applications
Example 1 :
Example 1: The Robertson Problem (Chemical Kinetics) This is a classic stiff DAE representing the reaction of three species. It is an Index-1 DAE where the total mass is conserved via an algebraic constraint.
double[] robertson_f(double t, double[] y) =>
[(-0.04 * y[0] + 1e4 * y[1] * y[2]),
(0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1]*y[1]),
y[0] + y[1] + y[2] - 1.0];
double[,] mass_f(double t, double[] y) => Diag([1, 1, 0]);
double[] y0 = [1.0, 0.0, 0.0];
(ColVec T, Matrix Y) = Ode45a(robertson_f, mass_f, y0, [0, 1e7]);
// Plot the result
Y[.., 1] = 1e4*Y[.., 1];
SemiLogx(T, Y);
Xlabel("Time t"); Ylabel("Soluton y");
Legend(["y_1", "1e4*y_2", "y_3"], MiddleLeft);
Title("Solution of Robertson's ODE with ODE45a");
SaveAs("Robertson-ODE-given-points-Ode45a.png");
Example 1 : The Simple Pendulum (Index-1)
A pendulum in Cartesian coordinates is naturally an Index-3 DAE. We solve the stabilized Index-1 version by including velocity constraints.
The position of the pendulum \((x, y)\) must satisfy the rigid rod constraint: \(x^2 + y^2 - 1 = 0\)
The Index-1 Formulation To reduce the index, we differentiate the constraint twice. The second derivative introduces the accelerations \(x''\) and \(y''\), allowing us to solve for the Lagrange multiplier \(\lambda\) (tension).
The resulting Index-1 system is:
\[\begin{split}\begin{array}{rcl}
x' &=& u \\
y' &=& v \\
u' &=& -\lambda x \\
v' &=& -\lambda y - g \\
0 &=& u^2 + v^2 - y g - \lambda
\end{array}\end{split}\]
double g = 9.81;
// State vector y = [x, y, u, v, λ]
double[] pendulum_f(double t, double[] y) =>
[y[2],
y[3],
-y[0] * y[4],
-y[1] * y[4] - g,
y[2]*y[2] + y[3]*y[3] - y[1] * g - y[4]];
double[,] mass_f(double t, double[] y) => Diag([1, 1, 1, 1, 0]);
double[] y0 = [0, 1, 1, 0, 1 - g];
var opts = Odeset(Stats: true, RelTol: 1e-6);
(ColVec T, Matrix Y) = Ode45a(pendulum_f, mass_f, y0, [0, 6], opts);
Plot(T, Y, Linewidth: 2); Xlabel("x"); Ylabel("y");
Legend(["x", "y", "u", "v", "λ"]);
Title("Pendulum Trajectory (DAE)");
SaveAs("Index_1-Pendulum-Problem-Ode45a.png");
Ouput
Summary of statistics by Ode45a
1321 successful steps
14 failed attempts
49672 function evaluations
5340 partial derivatives
5340 LU decompositions
17616 solutions of linear systems
As an exercise, the reader is encouraged to solve the problem using this initial condition y0 = [1, 0, 0, 1, 1];
Example 2 : Semi-Explicit DAE (The Transistor Amplifier)**
This example mimics the “hbdae” problem from MathWorks, representing an electrical circuit with nonlinear components.
The transistor amplifier circuit contains six resistors, three capacitors, and a transistor.
The initial voltage signal is \(U_e(t) = 0.4\sin(200\pi t)\).
The operating voltage is \(U_b = 6\).
The voltages at the nodes are given by: \(U_i(t)(i = 1, 2, 3, 4, 5)\).
The values of the resistors \(R_i(t)(i = 1, 2, 3, 4, 5)\). are constant, and the current through each resistor satisfies \(I = U/R\).
The values of the capacitors \(C_i(i = 1, 2, 3)\) are constant, and the current through each capacitor satisfies \(I=C⋅dU/dt\).
The goal is to solve for the output voltage through node 5, \(U_5(t)\). Using Kirchoff’s law to equalize the current through each node (1 through 5), you can obtain a system of five equations describing the circuit:
Node 1: \(C_1(U'_2 - U'_1) = (U_1 - U_e(t))/R_0\)
Node 2: \(C_1(U'_1 - U'_2) = (U_2 - U_b)/R_1 + U_2/R_1 + 0.01f(U_2 - U_3)\)
Node 3: \(-C_2U'_3 = U_3/R_3 - f(U_2 - U_3)\)
Node 4: \(C_3(U'_5 - U'_4) = (U_4 - U_b)/R_4 + 0.99f(U_2 - U_3)\)
Node 5: \(C_3(U'_4 - U'_5) = U_5/R_5\)
By extracting the coeeficients of the derivatives into a matrix, we have:
\[\begin{split}\begin{pmatrix}
-c_{1} & c_{1} & 0 & 0 & 0 \\
c_{1} & -c_{1} & 0 & 0 & 0 \\
0 & 0 & -c_{ 2} & 0 & 0 \\
0 & 0 & 0 & -c_{ 3} & c_{ 3} \\
0 & 0 & 0 & c_{ 3} & -c_{ 3}
\end{pmatrix}
\begin{pmatrix}
U'_1 \\ U'_2 \\ U'_3 \\ U'_4 \\ U'_5
\end{pmatrix} =
\begin{pmatrix}
(U_1 - U_e(t))/R_0 \\
(U_2 - U_b)/R_1 + U_2/R_1 + 0.01f(U_2 - U_3) \\
U_3/R_3 - f(U_2 - U_3) \\
(U_4 - U_b)/R_4 + 0.99f(U_2 - U_3) \\
U_5/R_5
\end{pmatrix}\end{split}\]
double Ub = 6, R0 = 1000, R15 = 9000, alpha = 0.99,
beta = 1e-6, Uf = 0.026, c1 = 1e-6, c2 = 2e-6, c3 = 3e-6;
double[,] Mass(double t, double[] y) => new double[,]
{
{-c1, c1, 0, 0, 0 },
{ c1, -c1, 0, 0, 0 },
{ 0, 0, -c2, 0, 0 },
{ 0, 0, 0, -c3, c3},
{ 0, 0, 0, c3, -c3}
};
double Ue(double t) => 0.4 * Sin(200 * pi * t);
double[] dudt(double t, double[] u)
{
double f23 = beta * (Exp((u[1] - u[2]) / Uf) - 1);
return [ -(Ue(t) - u[0])/R0,
-(Ub/R15 - u[1]*2/R15 - (1-alpha)*f23),
-(f23 - u[2]/R15),
-((Ub - u[3])/R15 - alpha*f23),
u[4]/R15 ];
}
double[] tspan = [0, 0.1];
double[] y0 = [0, Ub / 2, Ub / 2, Ub, 0];
var opts = Odeset(RelTol: 1e-5);
(ColVec T, Matrix Y) = Ode45a(dudt, Mass, y0, tspan, opts);
Scatter(T, Arrayfun(Ue, T), "o"); HoldOn();
Plot(T, Y[.., 4], "--r"); HoldOff();
Legend(["Input", "Output"], UpperLeft);
Xlabel("Time t"); Ylabel("Solution y");
Title("One Transistor Amplifier DAE Problem-Ode45a");
SaveAs("One-Transistor-Amplifier-DAE-Problem-Ode45a.png");
Example 3 : The Akzo Nobel Problem
A high-dimensional DAE describing a chemical process with 6 differential and 2 algebraic equations. This tests the solver’s ability to handle stiff systems with coupled variables.
Mathematical Description: The system is defined by reaction rates \(r_i\) and concentrations \(y_1, ..., y_8\):
\[\begin{split}\begin{array}{rcl}
r_1 &=& k_1 \cdot y_1^4 \cdot y_2^{0.5}\\
r_2 &=& k_2 \cdot y_3 \cdot y_4 \\
r_3 &=& k_2 / K \cdot y_1 \cdot y_5\\
r_4 &=& k_3 \cdot y_1 \cdot y_4^2\\
r_5 &=& k_4 \cdot y_6^2 \cdot y_2^{0.5}
\end{array}\end{split}\]
The differential equations are:
\[\begin{split}\begin{array}{rcl}
y_1' &=& -2r_1 + r_2 - r_3 - r_4\\
y_2' &=& -0.5r_1 - r_5 + 0.5F_{in}\\
y_3' &=& r_1 - r_2 + r_3\\
y_4' &=& -r_2 + r_3 - 2r_4\\
y_5' &=& r_2 - r_3 + r_4\\
y_6' &=& -r_5
\end{array}\end{split}\]
The algebraic constraints (Equilibrium):
\[\begin{split}\begin{array}{rcl}
0 &=& y_1 \cdot y_3 - y_7\\
0 &=& y_4 \cdot y_5 - y_8
\end{array}\end{split}\]
double k1 = 18.7, k2 = 0.58, k3 = 0.09, k4 = 0.42, K = 34.4, Fin = 0.012;
double r1(double[] y) => k1 * Pow(y[0], 4) * Pow(y[1], 0.5);
double r2(double[] y) => k2 * y[2] * y[3];
double r3(double[] y) => (k2 / K) * y[0] * y[4];
double r4(double[] y) => k3 * y[0] * Pow(y[3], 2);
double r5(double[] y) => k4 * Pow(y[5], 2) * Pow(y[1], 0.5);
double[] akzo_f(double t, double[] y) =>
[
-2*r1(y) + r2(y) - r3(y) - r4(y),
-0.5*r1(y) - r5(y) + 0.5*Fin,
r1(y) - r2(y) + r3(y),
-r2(y) + r3(y) - 2*r4(y),
r2(y) - r3(y) + r4(y),
-r5(y),
y[0] * y[2] - y[6],
y[3] * y[4] - y[7]
];
double[,] mass_f(double t, double[] y) => Diag([1, 1, 1, 1, 1, 1, 0, 0]);
double[] y0 = [0.444, 0.0012, 0.0, 0.0037, 0.0, 0.0, 0.0, 0.0];
(ColVec T, Matrix Y) = Ode45a(akzo_f, mass_f, y0, [0, 180]);
Plot(T, Y);
Xlabel("Time"); Ylabel("Concentration");
Title("Akzo Nobel Chemical Kinetics (DAE)");
SaveAs("Akzo-Nobel-Ode45a.png");
Index-2 DAE
Most DAE solvers usually avoid solving DAEs in index 2 form. But SepalSolver is able to handle most index 2 DAEs to a relative tolerance of \(10^{-4}\).
Now we look at examples of index 2 DAEs
Example 4 :
Usnig the example from “On the numerical solution of differential–algebraic equations with index-2” by Ercan Celık
\[\begin{split}\begin{align}
x'_1 &= \left(\alpha - \cfrac{1}{2 - t}\right)x_1 + (2 - t)\alpha z + \cfrac{3 - t}{2 - t}x_2 \\
x'_2 &= \cfrac{1 - \alpha}{t - 2} x_1 - x_2 + (\alpha - 1)z + 2e^t \\
0 &= (t + 2)x_1 + (t^2 - 4)x_2 - (t^2 + t - 2)e^t
\end{align}\end{split}\]
Intial condition: \(x_1(0) = 1, x_2(0) = 1\);
SepalSolver has the ability to compute consistent initial conditions for index 2 DAEs, so we can solve this problem without manually differentiating the algebraic constraint.
// define the DAE
double alpha = 10;
double[] Ercan(double t, double[] x) =>
[ (alpha - 1/(2-t))*x[0] + (2-t)*alpha*x[2] + (3-t)/(2-t)*x[1],
(1-alpha)/(t-2)*x[0] - x[1] + (alpha-1)*x[2] + 2*Exp(t),
(t+2)*x[0] + (t*t-4)*x[1] - (t*t+t-2)*Exp(t) ];
double[,] mass_f(double t, double[] x) => Diag([1, 1, 0]);
double[] y0 = [1, 1, 0]; // only the differential variables need initial conditions
var opts = Odeset(Stats: true);
(ColVec T, Matrix Y) = Ode45a(Ercan, mass_f, y0, [0, 1], opts);
Scatter(T, Hcart(Exp(T), Exp(T), -Exp(T).Div(2-T)), "o"); HoldOn();
Plot(T, Y); HoldOff();
Xlabel("Time t"); Ylabel("Solution x");
Legend(["x_1_Exact", "x_2_Exact", "z_Exact", "x_1_NumSol", "x_2_NumSol", "z_NumSol"]);
Title("Index-2 DAE Example (Ercan Celık)");
SaveAs("Index-2-DAE-Ercan-Celik.png");
// We can actually print out the result to compare with the analytical solution
Console.WriteLine("""
t || x_1_NumSol(t) | x_1_Exact(t) || x_2_NumSol(t) | x_2_Exact(t) || z_NumSol(t) | z_Exact(t)
--------++-----------------+----------------++-----------------+----------------++-----------------+---------------
""");
for (int i = 0; i < T.Numel; i++)
{
Console.WriteLine($"""
{T[i]:F2} || {Y[i, 0]:F6} | {Exp(T[i]):F6} || {Y[i, 1]:F6} | {Exp(T[i]):F6} || {Y[i, 2]:F6} | {-Exp(T[i])/(2-T[i]):F6}
""");
}
// We can compute the solution to a higher accuracy
Console.WriteLine("\n\nNow we compute the solution to a higher accuracy (RelTol = 1e-5):\n");
opts = Odeset(Stats: true, RelTol: 1e-5);
(T, Y) = Ode45a(Ercan, mass_f, y0, [0, 1], opts);
Console.WriteLine("""
t || x_1_NumSol(t) | x_1_Exact(t) || x_2_NumSol(t) | x_2_Exact(t) || z_NumSol(t) | z_Exact(t)
--------++-----------------+----------------++-----------------+----------------++-----------------+---------------
""");
for (int i = 0; i < T.Numel; i++)
{
Console.WriteLine($"""
{T[i]:F2} || {Y[i, 0]:F6} | {Exp(T[i]):F6} || {Y[i, 1]:F6} | {Exp(T[i]):F6} || {Y[i, 2]:F6} | {-Exp(T[i])/(2-T[i]):F6}
""");
}
Ouput
Summary of statistics by Ode45a
13 successful steps
0 failed attempts
355 function evaluations
52 partial derivatives
52 LU decompositions
135 solutions of linear systems
t || x_1_NumSol(t) | x_1_Exact(t) || x_2_NumSol(t) | x_2_Exact(t) || z_NumSol(t) | z_Exact(t)
--------++-----------------+----------------++-----------------+----------------++-----------------+---------------
0.00 || 1.000000 | 1.000000 || 1.000000 | 1.000000 || -0.500000 | -0.500000
0.01 || 1.010050 | 1.010050 || 1.010050 | 1.010050 || -0.507565 | -0.507563
0.04 || 1.038900 | 1.038900 || 1.038900 | 1.038900 || -0.529574 | -0.529555
0.09 || 1.096642 | 1.096647 || 1.096644 | 1.096647 || -0.574914 | -0.574840
0.18 || 1.193049 | 1.193078 || 1.193062 | 1.193078 || -0.654460 | -0.654292
0.28 || 1.318500 | 1.318555 || 1.318523 | 1.318555 || -0.765277 | -0.765061
0.38 || 1.457170 | 1.457228 || 1.457193 | 1.457228 || -0.897806 | -0.897604
0.48 || 1.610430 | 1.610486 || 1.610449 | 1.610486 || -1.057315 | -1.057121
0.58 || 1.779810 | 1.779863 || 1.779826 | 1.779863 || -1.250552 | -1.250374
0.68 || 1.967006 | 1.967052 || 1.967017 | 1.967052 || -1.486444 | -1.486291
0.78 || 2.173890 | 2.173929 || 2.173897 | 2.173929 || -1.776979 | -1.776864
0.88 || 2.402535 | 2.402563 || 2.402538 | 2.402563 || -2.138592 | -2.138532
0.98 || 2.655228 | 2.655243 || 2.655229 | 2.655243 || -2.594355 | -2.594369
1.00 || 2.718274 | 2.718282 || 2.718274 | 2.718282 || -2.718252 | -2.718282
Now we compute the solution to a higher accuracy (RelTol = 1e-5):
Summary of statistics by Ode45a
52 successful steps
0 failed attempts
1283 function evaluations
208 partial derivatives
208 LU decompositions
439 solutions of linear systems
t || x_1_NumSol(t) | x_1_Exact(t) || x_2_NumSol(t) | x_2_Exact(t) || z_NumSol(t) | z_Exact(t)
--------++-----------------+----------------++-----------------+----------------++-----------------+---------------
0.00 || 1.000000 | 1.000000 || 1.000000 | 1.000000 || -0.500000 | -0.500000
0.01 || 1.010050 | 1.010050 || 1.010050 | 1.010050 || -0.507565 | -0.507563
0.02 || 1.021438 | 1.021438 || 1.021438 | 1.021438 || -0.516197 | -0.516194
0.03 || 1.033831 | 1.033831 || 1.033831 | 1.033831 || -0.525663 | -0.525660
0.05 || 1.046984 | 1.046984 || 1.046984 | 1.046984 || -0.535796 | -0.535792
0.06 || 1.060740 | 1.060740 || 1.060740 | 1.060740 || -0.546487 | -0.546482
0.07 || 1.074997 | 1.074997 || 1.074997 | 1.074997 || -0.557668 | -0.557663
0.09 || 1.089691 | 1.089691 || 1.089691 | 1.089691 || -0.569300 | -0.569295
0.10 || 1.104787 | 1.104787 || 1.104787 | 1.104787 || -0.581365 | -0.581360
0.11 || 1.120271 | 1.120271 || 1.120271 | 1.120271 || -0.593863 | -0.593858
0.13 || 1.136135 | 1.136135 || 1.136135 | 1.136135 || -0.606796 | -0.606791
0.14 || 1.152381 | 1.152381 || 1.152381 | 1.152381 || -0.620175 | -0.620170
0.16 || 1.169015 | 1.169015 || 1.169015 | 1.169015 || -0.634017 | -0.634012
0.17 || 1.186048 | 1.186048 || 1.186048 | 1.186048 || -0.648341 | -0.648336
0.19 || 1.203493 | 1.203493 || 1.203493 | 1.203493 || -0.663171 | -0.663165
0.20 || 1.221361 | 1.221361 || 1.221361 | 1.221361 || -0.678527 | -0.678521
0.21 || 1.239669 | 1.239669 || 1.239669 | 1.239669 || -0.694438 | -0.694432
0.23 || 1.258435 | 1.258435 || 1.258435 | 1.258435 || -0.710934 | -0.710928
0.25 || 1.277679 | 1.277679 || 1.277679 | 1.277679 || -0.728047 | -0.728041
0.26 || 1.297424 | 1.297424 || 1.297424 | 1.297424 || -0.745815 | -0.745809
0.28 || 1.317691 | 1.317691 || 1.317691 | 1.317691 || -0.764276 | -0.764269
0.29 || 1.338503 | 1.338503 || 1.338503 | 1.338503 || -0.783468 | -0.783461
0.31 || 1.359885 | 1.359885 || 1.359885 | 1.359885 || -0.803437 | -0.803430
0.32 || 1.381867 | 1.381867 || 1.381867 | 1.381867 || -0.824232 | -0.824225
0.34 || 1.404478 | 1.404478 || 1.404478 | 1.404478 || -0.845908 | -0.845901
0.36 || 1.427755 | 1.427755 || 1.427755 | 1.427755 || -0.868526 | -0.868519
0.37 || 1.451730 | 1.451730 || 1.451730 | 1.451730 || -0.892148 | -0.892140
0.39 || 1.476441 | 1.476441 || 1.476441 | 1.476441 || -0.916844 | -0.916836
0.41 || 1.501933 | 1.501933 || 1.501933 | 1.501933 || -0.942695 | -0.942687
0.42 || 1.528252 | 1.528252 || 1.528252 | 1.528252 || -0.969788 | -0.969780
0.44 || 1.555453 | 1.555453 || 1.555453 | 1.555453 || -0.998224 | -0.998216
0.46 || 1.583591 | 1.583592 || 1.583591 | 1.583592 || -1.028111 | -1.028103
0.48 || 1.612734 | 1.612734 || 1.612734 | 1.612734 || -1.059576 | -1.059567
0.50 || 1.642953 | 1.642953 || 1.642953 | 1.642953 || -1.092758 | -1.092749
0.52 || 1.674327 | 1.674328 || 1.674327 | 1.674328 || -1.127815 | -1.127806
0.53 || 1.706948 | 1.706949 || 1.706949 | 1.706949 || -1.164930 | -1.164920
0.55 || 1.740913 | 1.740914 || 1.740913 | 1.740914 || -1.204303 | -1.204293
0.57 || 1.776352 | 1.776352 || 1.776352 | 1.776352 || -1.246191 | -1.246180
0.60 || 1.813401 | 1.813401 || 1.813401 | 1.813401 || -1.290875 | -1.290865
0.62 || 1.852219 | 1.852219 || 1.852219 | 1.852219 || -1.338691 | -1.338680
0.64 || 1.892995 | 1.892996 || 1.892995 | 1.892996 || -1.390040 | -1.390029
0.66 || 1.935966 | 1.935967 || 1.935967 | 1.935967 || -1.445418 | -1.445406
0.68 || 1.981419 | 1.981419 || 1.981419 | 1.981419 || -1.505437 | -1.505424
0.71 || 2.029695 | 2.029695 || 2.029695 | 2.029695 || -1.570846 | -1.570833
0.73 || 2.081236 | 2.081236 || 2.081236 | 2.081236 || -1.642613 | -1.642599
0.76 || 2.136636 | 2.136636 || 2.136636 | 2.136636 || -1.722042 | -1.722028
0.79 || 2.196676 | 2.196677 || 2.196677 | 2.196677 || -1.810879 | -1.810865
0.82 || 2.262492 | 2.262492 || 2.262492 | 2.262492 || -1.911658 | -1.911643
0.85 || 2.335823 | 2.335823 || 2.335823 | 2.335823 || -2.028282 | -2.028266
0.88 || 2.419602 | 2.419603 || 2.419603 | 2.419603 || -2.167350 | -2.167334
0.92 || 2.519828 | 2.519829 || 2.519828 | 2.519829 || -2.342279 | -2.342264
0.98 || 2.655280 | 2.655282 || 2.655281 | 2.655282 || -2.594455 | -2.594445
1.00 || 2.718281 | 2.718282 || 2.718281 | 2.718282 || -2.718279 | -2.718282
Example 5 : Pendulum position constraint (Index-2)
To reduce the index, if we differentiated the constraint once instead of twice, we end up with index 2 problem.
The resulting Index-1 system is:
\[\begin{split}\begin{array}{rcl}
x' &=& u \\
y' &=& v \\
u' &=& -\lambda x \\
v' &=& -\lambda y - g \\
0 &=& x u + y v
\end{array}\end{split}\]
double g = 9.81;
// State vector y = [x, y, u, v, λ]
double[] pendulum_f(double t, double[] y) =>
[y[2],
y[3],
-y[0] * y[4],
-y[1] * y[4] - g,
y[0]*y[2] + y[1]*y[3]];
double[,] mass_f(double t, double[] y) => Diag([1, 1, 1, 1, 0]);
double[] y0 = [0, 1, 1, 0, -1];
var opts = Odeset(Stats: true);
(ColVec T, Matrix Y) = Ode45a(pendulum_f, mass_f, y0, [0, 6], opts);
Plot(T, Y, Linewidth: 2); Xlabel("x"); Ylabel("y");
Legend(["x", "y", "u", "v", "λ"]);
Title("Pendulum Trajectory (DAE)");
SaveAs("Index_2-Pendulum-Problem-Ode45a.png");
Console.WriteLine("\n\n");
Console.WriteLine(Hcart(T, Y));
Ouput
Summary of statistics by Ode45a
240 successful steps
16 failed attempts
10496 function evaluations
1024 partial derivatives
1024 LU decompositions
4336 solutions of linear systems
0.0000 -0.0000 1.0000 1.0000 0.0000 -8.8100
0.0600 0.0603 0.9982 1.0159 -0.0614 -8.7513
0.1235 0.1263 0.9920 1.0671 -0.1358 -8.5678
0.1873 0.1969 0.9804 1.1534 -0.2316 -8.2253
0.2477 0.2699 0.9629 1.2658 -0.3548 -7.7075
0.3024 0.3423 0.9396 1.3890 -0.5060 -7.0207
0.3500 0.4113 0.9115 1.5078 -0.6804 -6.1938
0.3908 0.4749 0.8801 1.6116 -0.8697 -5.2688
0.4253 0.5320 0.8468 1.6950 -1.0649 -4.2904
0.4544 0.5822 0.8131 1.7567 -1.2579 -3.3000
0.4788 0.6256 0.7802 1.7983 -1.4420 -2.3339
0.4991 0.6624 0.7492 1.8230 -1.6117 -1.4238
0.5122 0.6864 0.7272 1.8328 -1.7300 -0.7796
0.5237 0.7075 0.7068 1.8367 -1.8385 -0.1783
0.5288 0.7169 0.6972 1.8368 -1.8885 0.1010
0.5330 0.7246 0.6892 1.8362 -1.9307 0.3385
0.5373 0.7325 0.6808 1.8347 -1.9739 0.5839
0.5421 0.7412 0.6713 1.8322 -2.0231 0.8649
0.5476 0.7513 0.6600 1.8280 -2.0809 1.1975
0.5540 0.7629 0.6465 1.8212 -2.1492 1.5944
0.5613 0.7763 0.6304 1.8109 -2.2297 2.0672
0.5697 0.7913 0.6114 1.7955 -2.3239 2.6277
0.5791 0.8082 0.5890 1.7732 -2.4333 3.2887
0.5897 0.8267 0.5626 1.7415 -2.5590 4.0640
0.6014 0.8468 0.5319 1.6974 -2.7026 4.9697
0.6142 0.8683 0.4961 1.6368 -2.8650 6.0238
0.6283 0.8908 0.4545 1.5548 -3.0472 7.2476
0.6436 0.9137 0.4063 1.4451 -3.2496 8.6656
0.6602 0.9365 0.3506 1.2997 -3.4717 10.3068
0.6781 0.9582 0.2861 1.1084 -3.7119 12.2052
0.6975 0.9774 0.2115 0.8584 -3.9666 14.4017
0.7186 0.9922 0.1251 0.5332 -4.2287 16.9455
0.7416 0.9997 0.0248 0.1113 -4.4857 19.8980
0.7670 0.9958 -0.0921 -0.4363 -4.7157 23.3389
0.7956 0.9733 -0.2295 -1.1505 -4.8786 27.3805
0.8291 0.9192 -0.3939 -2.0975 -4.8939 32.2129
0.8683 0.8137 -0.5813 -3.2896 -4.6048 37.7096
0.9114 0.6442 -0.7648 -4.5651 -3.8453 43.0818
0.9548 0.4219 -0.9066 -5.6188 -2.6147 47.2275
0.9938 0.1896 -0.9818 -6.2006 -1.1977 49.4396
1.0284 -0.0281 -0.9996 -6.3402 0.1785 49.9734
1.0596 -0.2238 -0.9746 -6.1439 1.4112 49.2509
1.0889 -0.3974 -0.9176 -5.7026 2.4700 47.5838
1.1171 -0.5503 -0.8349 -5.0785 3.3477 45.1584
1.1454 -0.6833 -0.7301 -4.3156 4.0393 42.0806
1.1745 -0.7962 -0.6049 -3.4480 4.5384 38.4045
1.2053 -0.8880 -0.4597 -2.5028 4.8346 34.1383
1.2396 -0.9563 -0.2921 -1.4993 4.9092 29.2135
1.2812 -0.9958 -0.0910 -0.4306 4.7135 23.3101
1.3278 -0.9929 0.1186 0.5074 4.2466 17.1606
1.3742 -0.9532 0.3021 1.1579 3.6541 11.7697
1.4137 -0.8999 0.4360 1.5146 3.1260 7.8183
1.4459 -0.8479 0.5301 1.6947 2.7107 5.0388
1.4722 -0.8022 0.5970 1.7818 2.3944 3.0632
1.4935 -0.7637 0.6455 1.8207 2.1541 1.6302
1.5071 -0.7389 0.6737 1.8330 2.0105 0.7954
1.5164 -0.7218 0.6921 1.8366 1.9155 0.2543
1.5224 -0.7107 0.7034 1.8369 1.8561 -0.0805
1.5274 -0.7015 0.7126 1.8362 1.8077 -0.3508
1.5322 -0.6927 0.7212 1.8346 1.7622 -0.6033
1.5375 -0.6831 0.7303 1.8319 1.7136 -0.8714
1.5435 -0.6722 0.7403 1.8278 1.6595 -1.1677
1.5503 -0.6596 0.7516 1.8217 1.5989 -1.4974
1.5583 -0.6452 0.7639 1.8130 1.5313 -1.8618
1.5673 -0.6289 0.7775 1.8011 1.4569 -2.2597
1.5776 -0.6104 0.7920 1.7853 1.3760 -2.6879
1.5892 -0.5899 0.8075 1.7650 1.2894 -3.1422
1.6022 -0.5671 0.8236 1.7396 1.1978 -3.6170
1.6166 -0.5422 0.8402 1.7086 1.1025 -4.1063
1.6327 -0.5150 0.8571 1.6716 1.0044 -4.6034
1.6504 -0.4857 0.8741 1.6285 0.9050 -5.1015
1.6700 -0.4543 0.8908 1.5792 0.8054 -5.5937
1.6916 -0.4208 0.9071 1.5241 0.7071 -6.0732
1.7154 -0.3853 0.9227 1.4637 0.6112 -6.5333
1.7416 -0.3478 0.9375 1.3989 0.5190 -6.9679
1.7705 -0.3084 0.9512 1.3308 0.4315 -7.3710
1.8024 -0.2670 0.9637 1.2614 0.3495 -7.7371
1.8379 -0.2236 0.9747 1.1927 0.2736 -8.0606
1.8773 -0.1779 0.9840 1.1279 0.2039 -8.3361
1.9212 -0.1296 0.9915 1.0709 0.1400 -8.5570
1.9702 -0.0783 0.9969 1.0269 0.0806 -8.7149
2.0247 -0.0231 0.9997 1.0028 0.0232 -8.7970
2.0844 0.0367 0.9993 1.0064 -0.0370 -8.7844
2.1457 0.0993 0.9950 1.0426 -0.1040 -8.6581
2.2090 0.1674 0.9859 1.1144 -0.1892 -8.3864
2.2707 0.2392 0.9710 1.2168 -0.2997 -7.9455
2.3277 0.3118 0.9501 1.3367 -0.4387 -7.3310
2.3781 0.3823 0.9240 1.4585 -0.6034 -6.5623
2.4217 0.4482 0.8939 1.5693 -0.7868 -5.6765
2.4587 0.5080 0.8613 1.6617 -0.9801 -4.7180
2.4900 0.5612 0.8277 1.7325 -1.1747 -3.7298
2.5163 0.6074 0.7944 1.7826 -1.3631 -2.7503
2.5383 0.6470 0.7625 1.8142 -1.5395 -1.8130
2.5565 0.6802 0.7330 1.8309 -1.6989 -0.9484
2.5677 0.7007 0.7135 1.8361 -1.8031 -0.3753
2.5748 0.7137 0.7004 1.8370 -1.8720 0.0086
2.5808 0.7247 0.6890 1.8363 -1.9315 0.3436
2.5860 0.7344 0.6787 1.8344 -1.9850 0.6478
2.5916 0.7446 0.6675 1.8310 -2.0424 0.9759
2.5977 0.7558 0.6548 1.8257 -2.1074 1.3517
2.6047 0.7685 0.6398 1.8173 -2.1829 1.7918
2.6126 0.7828 0.6222 1.8048 -2.2706 2.3099
2.6215 0.7989 0.6015 1.7863 -2.3724 2.9195
2.6315 0.8166 0.5772 1.7599 -2.4896 3.6347
2.6426 0.8359 0.5488 1.7228 -2.6239 4.4707
2.6549 0.8567 0.5158 1.6714 -2.7763 5.4447
2.6683 0.8787 0.4773 1.6014 -2.9481 6.5764
2.6829 0.9015 0.4328 1.5072 -3.1398 7.8887
2.6988 0.9245 0.3812 1.3818 -3.3516 9.4081
2.7161 0.9469 0.3215 1.2162 -3.5825 11.1659
2.7347 0.9676 0.2524 0.9991 -3.8302 13.1991
2.7549 0.9850 0.1725 0.7161 -4.0895 15.5523
2.7768 0.9968 0.0798 0.3484 -4.3513 18.2800
2.8009 0.9996 -0.0279 -0.1285 -4.5992 21.4510
2.8276 0.9881 -0.1539 -0.7483 -4.8042 25.1576
2.8582 0.9530 -0.3029 -1.5611 -4.9119 29.5392
2.8953 0.8752 -0.4838 -2.6550 -4.8026 34.8532
2.9367 0.7389 -0.6738 -3.9201 -4.2983 40.4224
2.9815 0.5349 -0.8449 -5.1530 -3.2625 45.4216
3.0232 0.3011 -0.9535 -5.9801 -1.8886 48.6063
3.0600 0.0733 -0.9973 -6.3221 -0.4645 49.9011
3.0928 -0.1341 -0.9909 -6.2721 0.8491 49.7274
3.1229 -0.3184 -0.9479 -5.9362 1.9941 48.4735
3.1516 -0.4811 -0.8766 -5.3908 2.9590 46.3839
3.1798 -0.6235 -0.7817 -4.6876 3.7392 43.5996
3.2083 -0.7460 -0.6659 -3.8647 4.3297 40.1974
3.2382 -0.8479 -0.5301 -2.9525 4.7224 36.2082
3.2706 -0.9276 -0.3735 -1.9745 4.9038 31.6069
3.3080 -0.9816 -0.1906 -0.9410 4.8449 26.2370
3.3520 -0.9998 0.0161 0.0726 4.5053 20.1695
3.3991 -0.9764 0.2157 0.8732 3.9529 14.3111
3.4417 -0.9281 0.3722 1.3584 3.3871 9.7014
3.4772 -0.8750 0.4840 1.6143 2.9187 6.4012
3.5061 -0.8262 0.5632 1.7424 2.5562 4.0606
3.5298 -0.7842 0.6204 1.8034 2.2795 2.3717
3.5491 -0.7491 0.6624 1.8291 2.0686 1.1327
3.5610 -0.7273 0.6862 1.8359 1.9460 0.4278
3.5686 -0.7134 0.7007 1.8371 1.8705 0.0009
3.5749 -0.7017 0.7124 1.8363 1.8087 -0.3445
3.5805 -0.6916 0.7222 1.8344 1.7564 -0.6349
3.5862 -0.6810 0.7322 1.8313 1.7032 -0.9284
3.5927 -0.6693 0.7430 1.8266 1.6454 -1.2441
3.5999 -0.6560 0.7547 1.8198 1.5818 -1.5895
3.6082 -0.6409 0.7675 1.8102 1.5116 -1.9671
3.6177 -0.6239 0.7814 1.7972 1.4349 -2.3760
3.6283 -0.6048 0.7963 1.7802 1.3521 -2.8133
3.6403 -0.5836 0.8120 1.7585 1.2639 -3.2747
3.6537 -0.5602 0.8283 1.7315 1.1711 -3.7547
3.6687 -0.5346 0.8450 1.6988 1.0748 -4.2471
3.6852 -0.5069 0.8620 1.6601 0.9762 -4.7453
3.7034 -0.4770 0.8789 1.6152 0.8766 -5.2425
3.7236 -0.4449 0.8955 1.5642 0.7772 -5.7318
3.7458 -0.4109 0.9116 1.5076 0.6795 -6.2064
3.7702 -0.3748 0.9270 1.4458 0.5845 -6.6599
3.7971 -0.3368 0.9415 1.3799 0.4935 -7.0861
3.8268 -0.2968 0.9549 1.3113 0.4075 -7.4792
3.8597 -0.2548 0.9669 1.2418 0.3272 -7.8336
3.8963 -0.2107 0.9775 1.1740 0.2531 -8.1441
3.9369 -0.1643 0.9863 1.1109 0.1851 -8.4047
3.9822 -0.1153 0.9933 1.0571 0.1227 -8.6086
4.0328 -0.0629 0.9980 1.0180 0.0642 -8.7463
4.0888 -0.0065 0.9999 1.0010 0.0065 -8.8036
4.1498 0.0548 0.9984 1.0139 -0.0557 -8.7592
4.2140 0.1213 0.9926 1.0629 -0.1299 -8.5846
4.2784 0.1923 0.9813 1.1477 -0.2249 -8.2509
4.3395 0.2656 0.9640 1.2594 -0.3470 -7.7408
4.3946 0.3384 0.9410 1.3828 -0.4974 -7.0603
4.4428 0.4079 0.9130 1.5025 -0.6712 -6.2375
4.4840 0.4719 0.8816 1.6073 -0.8603 -5.3145
4.5188 0.5294 0.8483 1.6918 -1.0558 -4.3363
4.5482 0.5800 0.8146 1.7545 -1.2492 -3.3445
4.5728 0.6237 0.7816 1.7971 -1.4340 -2.3758
4.5932 0.6608 0.7505 1.8224 -1.6046 -1.4620
4.6066 0.6851 0.7284 1.8327 -1.7239 -0.8127
4.6182 0.7064 0.7077 1.8368 -1.8335 -0.2055
4.6233 0.7159 0.6981 1.8371 -1.8838 0.0753
4.6277 0.7240 0.6898 1.8364 -1.9273 0.3205
4.6320 0.7319 0.6814 1.8350 -1.9710 0.5678
4.6368 0.7406 0.6719 1.8326 -2.0201 0.8485
4.6423 0.7507 0.6606 1.8284 -2.0776 1.1795
4.6486 0.7623 0.6472 1.8218 -2.1456 1.5738
4.6559 0.7755 0.6313 1.8117 -2.2256 2.0433
4.6642 0.7905 0.6124 1.7966 -2.3192 2.5999
4.6736 0.8073 0.5901 1.7746 -2.4278 3.2561
4.6841 0.8258 0.5640 1.7434 -2.5528 4.0260
4.6958 0.8458 0.5334 1.6999 -2.6955 4.9255
4.7086 0.8672 0.4979 1.6402 -2.8571 5.9725
4.7226 0.8896 0.4566 1.5593 -3.0384 7.1881
4.7378 0.9126 0.4087 1.4510 -3.2398 8.5967
4.7544 0.9354 0.3534 1.3074 -3.4610 10.2271
4.7722 0.9572 0.2893 1.1185 -3.7005 12.1130
4.7916 0.9765 0.2152 0.8716 -3.9546 14.2949
4.8126 0.9916 0.1294 0.5502 -4.2166 16.8217
4.8355 0.9995 0.0298 0.1333 -4.4743 19.7541
4.8608 0.9962 -0.0863 -0.4077 -4.7062 23.1707
4.8892 0.9749 -0.2227 -1.1131 -4.8735 27.1817
4.9224 0.9226 -0.3856 -2.0472 -4.8985 31.9715
4.9614 0.8199 -0.5725 -3.2312 -4.6274 37.4562
5.0043 0.6535 -0.7569 -4.5079 -3.8922 42.8540
5.0479 0.4321 -0.9017 -5.5819 -2.6749 47.0896
5.0873 0.1993 -0.9799 -6.1856 -1.2582 49.3868
5.1220 -0.0195 -0.9997 -6.3418 0.1240 49.9831
5.1534 -0.2163 -0.9762 -6.1572 1.3644 49.3044
5.1827 -0.3909 -0.9203 -5.7241 2.4310 47.6703
5.2110 -0.5446 -0.8386 -5.1062 3.3163 45.2719
5.2392 -0.6784 -0.7346 -4.3480 4.0156 42.2179
5.2682 -0.7921 -0.6102 -3.4838 4.5225 38.5643
5.2990 -0.8848 -0.4658 -2.5412 4.8269 34.3214
5.3331 -0.9541 -0.2992 -1.5398 4.9109 29.4253
5.3742 -0.9949 -0.0999 -0.4748 4.7280 23.5756
5.4206 -0.9939 0.1094 0.4701 4.2724 17.4351
5.4668 -0.9559 0.2934 1.1314 3.6855 12.0253
5.5064 -0.9032 0.4290 1.4988 3.1555 8.0262
5.5389 -0.8513 0.5246 1.6861 2.7359 5.2019
5.5653 -0.8053 0.5927 1.7775 2.4152 3.1905
5.5869 -0.7665 0.6422 1.8189 2.1710 1.7297
5.6007 -0.7412 0.6711 1.8324 2.0238 0.8725
5.6127 -0.7193 0.6946 1.8370 1.9023 0.1813
5.6180 -0.7095 0.7046 1.8371 1.8499 -0.1145
5.6224 -0.7015 0.7125 1.8364 1.8079 -0.3491
5.6268 -0.6934 0.7204 1.8349 1.7661 -0.5812
5.6317 -0.6843 0.7290 1.8325 1.7201 -0.8348
5.6374 -0.6739 0.7387 1.8287 1.6682 -1.1195
5.6441 -0.6618 0.7496 1.8230 1.6094 -1.4395
5.6518 -0.6478 0.7617 1.8149 1.5435 -1.7955
5.6606 -0.6319 0.7750 1.8037 1.4706 -2.1863
5.6706 -0.6138 0.7893 1.7887 1.3910 -2.6086
5.6819 -0.5937 0.8046 1.7693 1.3054 -3.0581
5.6947 -0.5713 0.8206 1.7448 1.2147 -3.5294
5.7088 -0.5468 0.8372 1.7149 1.1200 -4.0163
5.7246 -0.5200 0.8541 1.6790 1.0224 -4.5124
5.7420 -0.4911 0.8710 1.6371 0.9231 -5.0108
5.7613 -0.4601 0.8878 1.5889 0.8235 -5.5046
5.7825 -0.4270 0.9042 1.5349 0.7248 -5.9868
5.8058 -0.3918 0.9199 1.4754 0.6284 -6.4509
5.8315 -0.3547 0.9349 1.4113 0.5355 -6.8906
5.8599 -0.3156 0.9488 1.3438 0.4470 -7.2998
5.8913 -0.2746 0.9615 1.2745 0.3640 -7.6730
5.9260 -0.2315 0.9728 1.2055 0.2869 -8.0046
5.9647 -0.1862 0.9824 1.1398 0.2161 -8.2892
6.0000 -0.1469 0.9891 1.0905 0.1619 -8.4861
Observe that the initial condition supplied for \(\lambda\) was \(-1\); but the result returned shown that the correct initial condition for the algebraic variable \(\lambda\) is \(-8.81\). Sending in a wrong initial condition was done on purpose, to test the ability of sepalsolver to compute the initial condition of the algebraic variable.