dydx

Chapter 08

Ordinary Differential Equations

Your Progress in this Chapter

0%

0 / 21 units completed

Chapter Overview

Differential Algebraic Equations

Section 8.6 of 726 min7 code examples

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

.. Admonition:: 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.

Example 1C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Code is ready to run
OutputFrom the book
Robertson-ODE-given-points-Ode45a.png

.. Admonition:: 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:

x′=uy′=vu′=−λxv′=−λy−g0=u2+v2−yg−λ\begin{array}{rcl} x' &=& u \\ y' &=& v \\ u' &=& -\lambda x \\ v' &=& -\lambda y - g \\ 0 &=& u^2 + v^2 - y g - \lambda \end{array}

Example 2C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Code is ready to run
OutputFrom the book
Summary of statistics by Ode45a
        1320 successful steps
        19 failed attempts
        36985 function evaluations
        1339 partial derivatives
        5356 LU decompositions
        23577 solutions of linear systems
Index_1-Pendulum-Problem-Ode45a.png

As an exercise, the reader is encouraged to solve the problem using this initial condition y0 = [1, 0, 0, 1, 1];

.. Admonition:: 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:

(−c1c1000c1−c100000−c200000−c3c3000c3−c3)(U1′U2′U3′U4′U5′)=((U1−Ue(t))/R0(U2−Ub)/R1+U2/R1+0.01f(U2−U3)U3/R3−f(U2−U3)(U4−Ub)/R4+0.99f(U2−U3)U5/R5)\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}

Example 3C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
Code is ready to run
OutputFrom the book
One-Transistor-Amplifier-DAE-Problem-Ode45a.png

.. Admonition:: 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:

r1=k1⋅y14⋅y20.5r2=k2⋅y3⋅y4r3=k2/K⋅y1⋅y5r4=k3⋅y1⋅y42r5=k4⋅y62⋅y20.5\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}

The differential equations are:

y1′=−2r1+r2−r3−r4y2′=−0.5r1−r5+0.5Finy3′=r1−r2+r3y4′=−r2+r3−2r4y5′=r2−r3+r4y6′=−r5\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}

The algebraic constraints (Equilibrium):

0=y1⋅y3−y70=y4⋅y5−y8\begin{array}{rcl} 0 &=& y_1 \cdot y_3 - y_7\\ 0 &=& y_4 \cdot y_5 - y_8 \end{array}

Example 4C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
Code is ready to run
OutputFrom the book
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

.. Admonition:: Example 4 :

Usnig the example from "On the numerical solution of differential–algebraic equations with index-2" by Ercan Celık

x1′=(α−12−t)x1+(2−t)αz+3−t2−tx2x2′=1−αt−2x1−x2+(α−1)z+2et0=(t+2)x1+(t2−4)x2−(t2+t−2)et\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}

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.

Example 5C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
Code is ready to run
OutputFrom the book
Summary of statistics by Ode45a
        13 successful steps
        1 failed attempts
        581 function evaluations
        14 partial derivatives
        55 LU decompositions
        456 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.07  ||     1.068361    |    1.068367    ||     1.068364    |    1.068367    ||     -0.552529   |  -0.552451
  0.17  ||     1.180678    |    1.180728    ||     1.180701    |    1.180728    ||     -0.644069   |  -0.643846
  0.27  ||     1.304848    |    1.304907    ||     1.304873    |    1.304907    ||     -0.752803   |  -0.752599
  0.37  ||     1.442086    |    1.442145    ||     1.442109    |    1.442145    ||     -0.882858   |  -0.882657
  0.47  ||     1.593760    |    1.593817    ||     1.593780    |    1.593817    ||     -1.039277   |  -1.039083
  0.57  ||     1.761387    |    1.761440    ||     1.761403    |    1.761440    ||     -1.228633   |  -1.228453
  0.67  ||     1.946645    |    1.946692    ||     1.946657    |    1.946692    ||     -1.459589   |  -1.459433
  0.77  ||     2.151388    |    2.151428    ||     2.151395    |    2.151428    ||     -1.743763   |  -1.743644
  0.87  ||     2.377665    |    2.377695    ||     2.377669    |    2.377695    ||     -2.097043   |  -2.096976
  0.88  ||     2.401569    |    2.401591    ||     2.401571    |    2.401591    ||     -2.136823   |  -2.136897
  0.95  ||     2.577929    |    2.577937    ||     2.577929    |    2.577937    ||     -2.448177   |  -2.448159
  1.00  ||     2.718279    |    2.718282    ||     2.718279    |    2.718282    ||     -2.718276   |  -2.718282


Now we compute the solution to a higher accuracy (RelTol = 1e-5):

Summary of statistics by Ode45a
        28 successful steps
        3 failed attempts
        1238 function evaluations
        31 partial derivatives
        121 LU decompositions
        979 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.03  ||     1.032886    |    1.032886    ||     1.032886    |    1.032886    ||     -0.524948   |  -0.524936
  0.06  ||     1.065689    |    1.065689    ||     1.065689    |    1.065689    ||     -0.550377   |  -0.550352
  0.10  ||     1.104920    |    1.104922    ||     1.104921    |    1.104922    ||     -0.581502   |  -0.581469
  0.14  ||     1.148497    |    1.148500    ||     1.148498    |    1.148500    ||     -0.616999   |  -0.616961
  0.18  ||     1.195493    |    1.195496    ||     1.195495    |    1.195496    ||     -0.656388   |  -0.656347
  0.22  ||     1.245612    |    1.245615    ||     1.245613    |    1.245615    ||     -0.699681   |  -0.699638
  0.26  ||     1.298867    |    1.298871    ||     1.298869    |    1.298871    ||     -0.747165   |  -0.747120
  0.30  ||     1.355433    |    1.355437    ||     1.355435    |    1.355437    ||     -0.799302   |  -0.799255
  0.35  ||     1.415576    |    1.415580    ||     1.415578    |    1.415580    ||     -0.856699   |  -0.856650
  0.39  ||     1.479636    |    1.479640    ||     1.479637    |    1.479640    ||     -0.920110   |  -0.920059
  0.44  ||     1.548019    |    1.548024    ||     1.548021    |    1.548024    ||     -0.990458   |  -0.990405
  0.48  ||     1.621215    |    1.621219    ||     1.621216    |    1.621219    ||     -1.068881   |  -1.068827
  0.53  ||     1.699808    |    1.699813    ||     1.699810    |    1.699813    ||     -1.156799   |  -1.156743
  0.58  ||     1.784512    |    1.784517    ||     1.784513    |    1.784517    ||     -1.256006   |  -1.255948
  0.63  ||     1.876211    |    1.876216    ||     1.876212    |    1.876216    ||     -1.368817   |  -1.368758
  0.68  ||     1.976027    |    1.976033    ||     1.976029    |    1.976033    ||     -1.498292   |  -1.498233
  0.73  ||     2.085432    |    2.085438    ||     2.085433    |    2.085438    ||     -1.648598   |  -1.648540
  0.79  ||     2.206429    |    2.206434    ||     2.206430    |    2.206434    ||     -1.825633   |  -1.825578
  0.80  ||     2.219616    |    2.219621    ||     2.219617    |    2.219621    ||     -1.845573   |  -1.845588
  0.82  ||     2.269376    |    2.269379    ||     2.269376    |    2.269379    ||     -1.922397   |  -1.922397
  0.86  ||     2.364151    |    2.364153    ||     2.364151    |    2.364153    ||     -2.074605   |  -2.074582
  0.87  ||     2.377149    |    2.377151    ||     2.377149    |    2.377151    ||     -2.096067   |  -2.096072
  0.89  ||     2.432778    |    2.432779    ||     2.432778    |    2.432779    ||     -2.189791   |  -2.189788
  0.93  ||     2.542586    |    2.542587    ||     2.542586    |    2.542587    ||     -2.383354   |  -2.383338
  0.94  ||     2.558171    |    2.558172    ||     2.558171    |    2.558172    ||     -2.411757   |  -2.411761
  0.97  ||     2.625979    |    2.625980    ||     2.625979    |    2.625980    ||     -2.538294   |  -2.538292
  1.00  ||     2.718281    |    2.718282    ||     2.718281    |    2.718282    ||     -2.718283   |  -2.718282
Index-2-DAE-Ercan-Celik.png

.. Admonition:: 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:

x′=uy′=vu′=−λxv′=−λy−g0=xu+yv\begin{array}{rcl} x' &=& u \\ y' &=& v \\ u' &=& -\lambda x \\ v' &=& -\lambda y - g \\ 0 &=& x u + y v \end{array}

Example 6C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Code is ready to run
OutputFrom the book
Summary of statistics by Ode45a
        316 successful steps
        134 failed attempts
        15715 function evaluations
        450 partial derivatives
        1610 LU decompositions
        11388 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.1275    0.1306    0.9914    1.0715   -0.1411   -8.5506
   0.2089    0.2222    0.9750    1.1904   -0.2713   -8.0594
   0.3049    0.3459    0.9383    1.3951   -0.5142   -6.9615
   0.3164    0.3621    0.9322    1.4232   -0.5528   -6.8125
   0.3688    0.4401    0.8980    1.5559   -0.7626   -5.7912
   0.3780    0.4545    0.8908    1.5794   -0.8059   -5.5935
   0.4213    0.5253    0.8510    1.6858   -1.0405   -4.4095
   0.4289    0.5382    0.8429    1.7032   -1.0875   -4.1842
   0.4655    0.6019    0.7986    1.7771   -1.3395   -2.8690
   0.4721    0.6137    0.7896    1.7882   -1.3899   -2.6159
   0.5045    0.6722    0.7404    1.8277   -1.6595   -1.1564
   0.5104    0.6831    0.7304    1.8319   -1.7134   -0.8727
   0.5399    0.7373    0.6757    1.8335   -2.0007    0.7477
   0.5454    0.7473    0.6646    1.8299   -2.0577    1.0633
   0.5726    0.7966    0.6046    1.7892   -2.3575    2.8392
   0.5776    0.8056    0.5926    1.7771   -2.4160    3.1822
   0.6028    0.8493    0.5281    1.6914   -2.7202    5.0899
   0.6074    0.8571    0.5152    1.6704   -2.7790    5.4589
   0.6310    0.8949    0.4463    1.5372   -3.0825    7.4957
   0.6354    0.9017    0.4326    1.5067   -3.1407    7.8904
   0.6576    0.9333    0.3593    1.3236   -3.4379   10.0544
   0.6618    0.9387    0.3448    1.2833   -3.4943   10.4744
   0.6831    0.9637    0.2673    1.0481   -3.7786   12.7616
   0.6872    0.9678    0.2519    0.9975   -3.8319   13.2059
   0.7077    0.9854    0.1705    0.7086   -4.0956   15.6096
   0.7116    0.9881    0.1544    0.6474   -4.1441   16.0768
   0.7316    0.9976    0.0692    0.3039   -4.3784   18.5880
   0.7351    0.9986    0.0539    0.2383   -4.4167   19.0337
   0.7496    1.0000   -0.0114   -0.0523   -4.5654   20.9590
   0.7728    0.9929   -0.1197   -0.5736   -4.7582   24.1473
   0.7759    0.9910   -0.1341   -0.6468   -4.7786   24.5668
   0.7894    0.9800   -0.1993   -0.9871   -4.8536   26.4858
   0.8151    0.9457   -0.3250   -1.6889   -4.9139   30.1861
   0.8186    0.9397   -0.3421   -1.7884   -4.9127   30.6860
   0.8341    0.9084   -0.4183   -2.2456   -4.8770   32.9270
   0.8370    0.9018   -0.4323   -2.3321   -4.8646   33.3410
   0.8517    0.8642   -0.5032   -2.7784   -4.7722   35.4242
   0.8545    0.8562   -0.5168   -2.8660   -4.7483   35.8266
   0.8694    0.8102   -0.5863   -3.3227   -4.5917   37.8682
   0.8723    0.8004   -0.5996   -3.4120   -4.5545   38.2636
   0.8875    0.7449   -0.6672   -3.8739   -4.3249   40.2490
   0.8905    0.7332   -0.6801   -3.9633   -4.2728   40.6322
   0.9062    0.6676   -0.7446   -4.4195   -3.9621   42.5246
   0.9093    0.6538   -0.7567   -4.5066   -3.8934   42.8881
   0.9255    0.5771   -0.8167   -4.9437   -3.4936   44.6434
   0.9287    0.5611   -0.8278   -5.0256   -3.4067   44.9784
   0.9456    0.4725   -0.8814   -5.4270   -2.9090   46.5448
   0.9490    0.4539   -0.8911   -5.5003   -2.8020   46.8409
   0.9661    0.3570   -0.9342   -5.8298   -2.2278   48.0954
   0.9693    0.3385   -0.9410   -5.8826   -2.1163   48.3101
   0.9844    0.2477   -0.9689   -6.0991   -1.5592   49.1193
   0.9872    0.2308   -0.9730   -6.1317   -1.4544   49.2532
   1.0003    0.1497   -0.9888   -6.2549   -0.9473   49.7070
   1.0241   -0.0005   -1.0000   -6.3437    0.0034   50.0166
   1.0268   -0.0180   -0.9999   -6.3424    0.1141   50.0429
   1.0344   -0.0661   -0.9979   -6.3264    0.4191   49.9805
   1.0492   -0.1593   -0.9873   -6.2432    1.0071   49.6605
   1.0724   -0.3013   -0.9536   -5.9801    1.8892   48.6528
   1.0756   -0.3204   -0.9473   -5.9316    2.0060   48.4965
   1.0905   -0.4066   -0.9136   -5.6723    2.5244   47.4952
   1.0934   -0.4232   -0.9061   -5.6146    2.6220   47.2835
   1.1090   -0.5080   -0.8614   -5.2761    3.1119   45.9572
   1.1122   -0.5250   -0.8511   -5.1993    3.2074   45.6655
   1.1297   -0.6120   -0.7909   -4.7544    3.6789   43.8828
   1.1331   -0.6279   -0.7783   -4.6629    3.7616   43.5239
   1.1502   -0.7036   -0.7106   -4.1777    4.1365   41.5227

Showing the first 80 of 329 lines.

Index_2-Pendulum-Problem-Ode45a.png

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.

<header> Solving Index 3 </header> To show more capability of the sepalsolver with higher index DAEs, we present this solution of the Pendulum equation from index 0 to index 3 below

% --- Index 0 ---

x˙=uy˙=vu˙=−xλv˙=−yλ−gλ˙=−2λ(xu+yv)−3gv\begin{array}{rcl} \dot{x} &= u \\ \dot{y} &= v \\ \dot{u} &= -x \lambda \\ \dot{v} &= -y \lambda - g \\ \dot{\lambda} &= -2\lambda(xu + yv) - 3gv \end{array}

% --- Index 1 ---

x˙=uy˙=vu˙=−xλv˙=−yλ−g0=u2+v2−yg−λ\begin{array}{rcl} \dot{x} &= u \\ \dot{y} &= v \\ \dot{u} &= -x \lambda \\ \dot{v} &= -y \lambda - g \\ 0 &= u^2 + v^2 - y g - \lambda \end{array}

% --- Index 2 ---

x˙=uy˙=vu˙=−xλv˙=−yλ−g0=xu+yv\begin{array}{rcl} \dot{x} &= u \\ \dot{y} &= v \\ \dot{u} &= -x \lambda \\ \dot{v} &= -y \lambda - g \\ 0 &= x u + y v \end{array}

% --- Index 3 ---

x˙=uy˙=vu˙=−xλv˙=−yλ−g0=x2+y2−1\begin{array}{rcl} \dot{x} &= u \\ \dot{y} &= v \\ \dot{u} &= -x \lambda \\ \dot{v} &= -y \lambda - g \\ 0 &= x^2 + y^2 - 1 \end{array}

The result is assessed using these errors

r=∣x2+y2−1∣ϵ=∣xu+yv∣\begin{array}{rcl} r &= |x^2 + y^2 - 1| \\ \epsilon &= |xu + yv| \end{array}

Example 7C#

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Code is ready to run
OutputFrom the book
Summary of statistics by Ode45
        192 successful steps
        3 failed attempts
        1171 function evaluations

Summary of statistics by Ode45a
        1222 successful steps
        10 failed attempts
        34363 function evaluations
        1232 partial derivatives
        4928 LU decompositions
        22026 solutions of linear systems

Summary of statistics by Ode45a
        965 successful steps
        13 failed attempts
        35795 function evaluations
        978 partial derivatives
        3912 LU decompositions
        25998 solutions of linear systems

Summary of statistics by Ode45a
        10949 successful steps
        26 failed attempts
        343161 function evaluations
        10975 partial derivatives
        43897 LU decompositions
        233402 solutions of linear systems
Pendulum-Problem-Ode45a.png