146
Xnumbers Tutorial
223
ODE Implicit Predictor-Corrector
ODE_PC2I(Equations, VarInit, Step, [Par, …])
This function integrates numerically an ordinary differential equation of 1st order or a differential
system of 1st order, with the 2th order implicit Predictor-Corrector method. Thanks to its large
stability it is suitable for solving "stiff" problems.
The input arguments are identical to the function ODE_RK4
The function can return all the solution points in an array. Before inserting the function, select
as many rows that you want to fill. For example if you select a (p x n+1) array , the function will
return p solution points of a n x n differential system.
This function uses the 1
st
order Euler formula as predictor and the 2
nd
order trapezoidal formula
as the corrector.
( , , )
1
n
n
n
n
h f f t t y
y
y
+ ⋅
+
=
+
Predictor (Euler formula)
(
)
)
( ,
( , , )
/2
1
1
1
+
+
+
+
⋅
+
=
n
n
n
n
n
n
y
f t
f t t y
h
y
y
Corrector (Trapezoidal formula)
The second non-linear equation is started with the predicted value and then solved respect to
the variable yn+1 till the convergence.
This algorithm, despite its low order, exibits a large stability and, thus, it is suitable for solving
"stiff" problems
Let’s see how it works with an example
Solve numerically the following 2
nd
order Cauchy’s problem for 1 ≤ x ≤ 7 and k = 100
0
( 1)
2
2
=
+
+ +
ky
dx
dy
k
dx
d y
,
0.5
(1) 1 , '(1)
=−
=
y
y
First of all, we transform the problem in a 1st order differential equations system taking
'
,
2
1
y
y y
y
=
=
+
−
=−
=
2
1
2
2
1
( 1)
'
'
y
k
ky
y
y
y
=−
=
0.5
(1)
(1) 1
2
1
y
y
A possible arrangement may be the following
We have written in cell A2 and A3 the
differential equations
In cell D3 we have inserted the k parameter
and in the cell E3 the h step
In the range A8:C8 we have put the starting
values for x and y1 and y2.
note that we have also added the labels "k" ,
"x", "y1", "y2" just above theirs values.
Labels
are necessary for the correct
assignment in the symbolic equation.
Now select the range A9:C40 and insert the
function ODE_PC2I passing its arguments
The scatter plot of the 30 points (x, y1) is shown in the following graph
VB.NET PDF - How to Add Barcode on PDF Page text in PDF, C#.NET edit PDF bookmark, C#.NET edit PDF metadata, C#.NET VB.NET PDF barcode creator add-on, which combines the PDF reading add-on with
preview edit pdf metadata; google search pdf metadata
83
Xnumbers Tutorial
224
Compare the approximate solution (dotted line) with the exact one (pink continue line)
The fitting, even with only 30 points, looks good. If we try to solve this problem with the 4
th
order
Runge-Kutta algorithm we have to choose a step less then 0.025 for avoiding the instability
and, thus, we need more than 250 points for reaching the same integration interval.
The Lotcka-Volterra Model
The following two-dimensional differential system
+
=−
−
=
a y y a a xy
dt
dy
ax a a xy
dt
dx
4
3
2
1
is called the Lockta-Volterra model or also "prey-predator" model. It is very useful in biology,
chemistry and many other fields. The numerical integration of this ODE family requires a stable
algorithm otherwise the result are not acceptable
The function x(t) and y(t), depending from the time, may represent different things. In a
biological models, for example, they simulate respectively the population of pray (rabbits) and
predator (foxes) at time t. The proportionality constants a
1
, a
2
, a
3
, and a
4
are positive. It is
known that its solution leads to oscillate steady-state.
. Let's see with a practical example
=− +
−
=
xy
x
y
xy
x
x
2
' 5
2
' 2
=
=
(0) 1
(0) 1
y
x
Insert the equation defintion in the cells
A1 and A2 , the step in the cell C2 = 0.02
and the starting values [0, 1, 1] in the
range A6:C6
Do not miss the labels "t", "x", "y"
Then, select the range A7:C407 and
insert the function ODE_PC2I with the
ctrl-shift-enter keys sequence
The following graph at the left shows the result of the function x(t) and y(t). By performing the
scatter plot of x-y variable we have the steady-state plot at the right. It shows a closed loop, a
limit cycle.
VB.NET PDF - Convert CSV to PDF pages, C#.NET search text in PDF, C#.NET edit PDF bookmark, C#.NET edit PDF metadata, C#.NET VB.NET Demo Code for Converting RTF to PDF. Add necessary references
batch pdf metadata; modify pdf metadata
109
Xnumbers Tutorial
225
0
1
2
3
4
5
6
0
2
4
6
8
10
12
x
y
0
1
2
3
0
2
4
6
We can easily change the zoom of the graph by changing the integration step h
Note that the oscillations are quite deep. The above system has two equilibrium points where
dx/dt = 0 and dy/dt = 0 simultaneously, that can be found solving the algebric system
=
+
−
=
−
0
2
5
0
2
2
xy
x
xy
x
⇒
=
=
0
0
y
x
=
=
1
2.5
y
x
The non trivial solutions means that if the initial populations start from x
0
= 2.5 and y
0
= 1, then
there will be no oscillation and the population fraction x/y will be always constant in the time.
On the contrary, the populations become progressively to oscillate when the initial populations
are different from the equilibrium point [2.5 , 1].
Here same examples:
x(0) = 2.45 and y(0) = 1
x(0) = 2.4 and y(0) = 1
0
0.5
1
1.5
2
2.5
3
0
2
4
6
8
10
12
x
y
0
0.5
1
1.5
2
2.5
3
0
2
4
6
8
10
12
x
y
190
Xnumbers Tutorial
226
Differential Systems
OD Linear System
= ODE_SYSL(A, y0, step, [b]))
This function integrates numerically a system of ordinary linear differential equations of 1st
order with constant coefficients, starting from an initial value. For example, the general form of
a 3x3 linear differential system is:
+
+
+
=
+
+
+
=
+
+
+
=
3
1
33
1
32
1
31
3
2
1
23
1
22
1
21
2
1
1
13
1
12
1
11
1
'
'
'
a y y a a y y a a y y b
y
a y y a a y y a a y y b
y
a y y a a y y a a y y b
y
where all the coefficients aij and bi are constant.
Such system can be put in the following handy matrix form.
[ ]
]
A y y b
y
⋅ +
'=
where
=
33
32
31
23
22
21
13
12
11
a
a
a
a
a
a
a
a
a
A
[
]
T
y y y y
y
3
2
1
, ,
=
[
]
T
b b b b
b
3
2
1
= , ,
and with the initial condition:
[
]
T
y
y
y
y
(0)
3
(0)
2
(0)
1
0
,
,
=
The constant term b is optional. If omitted the system is called "homogeneous"
This function uses the exponential expansion method that, for this kind of differential systems,
is both accurate and efficient.
The function returns an (n x m ) array containing all the nodes of the integration: m is the
number of equations; n is the numbers of the nodes. The function automatically sets n equal to
the rows of the range that you have selected before inserting it.
Let's see how it works practically
Solve the following homogeneous differential system with constant coefficients
[ ]
]
A y
y
⋅
'=
'
where
−
−
−
−
−
=
1
2
2
16
6
16
19
10
20
A
[
]
T
y
1 , 0 , 0
0
=
For 0 ≤ x ≤ 4 and h = 0.05
The numerical solution can be arranged as in the following worksheet
The initial values are in the first row (range B7:D7). The column "x" was added only for clarity
but it is not indispensable at all. Select the range B8:D87, than insert the function ODE_SYSL
giving the suitable parameters A, y
0
, h. Then press CTRL+SHIFT+ENTER
315
Xnumbers Tutorial
227
All the 240 cells will be filled with the nodal solutions of y1, y2, y3.
The scale can be esily arranged simply by changing the parameter h
-1
-0.5
0
0.5
1
1.5
0
0.5
1
1.5
2
2.5
3
3.5
4
y1
y2
y3
-1
-0.5
0
0.5
1
1.5
0
0.25
0.5
0.75
1
y1
y2
y3
h = 0.05
h = 0.012
Compare with the exact solutions
+
=−
−
=
+
+
=−
−
−
−
−
−
−
−
x
x
x
x
x
x
x
e
e
y
e
e
y
e
e
e
y
2
3
10
2
2
10
2
1
2
2
2
2
2
2
High order linear ODE
The function ODE_SYSL can also be used to solve high order linear ODE with constant
coefficients, that in general can be written as
a y y ay y a a y y b
a y
y
n
n
n
=
+
+
+
+
−
−
0
1
2
( 1)
1
( )
'
''
....
with the initial conditions
0
0
( )
y
y x
=
,
'
'( )
0
0
y
y x
=
,
''
''( )
0
0
y
y x
=
, ,
( 1)
0
0
( 1)
( )
−
−
=
n
n
y
x
y
As known, such ODE can be transformed into a linear differential system of 1
st
order.
+
−
−
−
−
=
−
−
−
b
y
y
y
y
a
a
a
a
y
y
y
y
n
n
n
n
n
0
0
0
1
0
0
0
0
1
0
0
0
0
1
0
'
'
'
'
1
2
1
1
2
1
0
1
2
1
L
L
M
M
K
O
K
K
K
M
M
L
having the initial conditions
[
]
[
]
T
n
T
n
y
y
y
y
y
y
y
y
y
(0)
(0) , '(0) , ''(0) ,
(0) , (0) , (0) , (0)
1
3
2
1
0
−
=
=
219
Xnumbers Tutorial
228
For example assume tha you have to calculate the solution of following IVP problem:
0
''' 5 5 '' ' 9 9 ' ' 5
=
+
+
+
y
y
y
y
, with
(0 )=2
y
,
3
'(0)
=−
=
y
,
''(0) 4
=
y
Introducing the axiliary variables y
1
= y , y
2
= y' , y
3
= y'' , we get the following equivalent
differential system
⋅
−
−
−
=
3
2
1
3
2
1
5
9
5
1
0
0
0
1
0
'
'
'
y
y
y
y
y
y
with
=
=−
=
(0) 4
3
(0)
(0) 2
3
2
1
y
y
y
Observe that the last row of the matrix contains the coefficients of the given ODE with the
opposite sign; besides of that, it has only all "1" in the upper subdiagonal.
Insert the initial values in the first row (range B10:D10). The column "x" was added only for
clarity but it is not indispensable at all.Select the range B11:D40, where you want to otput the
results and insert the function ODE_SYSL giving the suitable parameters: A, y
0
, h.
Then press CTRL+SHIFT+ENTER
The selected area will be filled with the numerical solution of the given system
0
0.5
1
1.5
2
2.5
0
1
2
3
4
5
6
y1
Exact
Observe that the solution y
1
(x) is also the solution y(x) of the given ODE. Compare with the
exact solution
x
e
e
y x
x
x
cos
( )
−2
−
+
=
Of course the above differential system can be also solved with other methods.
In order to show the accuracy of the
exponential method we put in a graph the
average relative error
e(x) = |yi - y(xi)| / | y(xi)| ,
obtained in the same condition, with 3
different methods: Exponential, Euler,
and Runge-Kutta 4.
The graph is eloquent. The error of the
Exponetial method is several times more
accurate then the others
1E-16
1E-14
1E-12
1E-10
1E-08
1E-06
0.0001
0.01
1
100
0
1
2
3
4
5
6
Euler
RK 4
Expo
Clearly it takes advantages using dedicated methods for linear differential equations.
Another important feature of the exponential method is its high stability
Let' s try the following test stiff system
−
=−
=
2
1
2
2
1
21
20
'
y
y
y
y
y
=−
=
6770293
0.
(0.4)
6706555
(0.4) 0.
2
1
y
y
95
Xnumbers Tutorial
229
The exact solution is
−
=−
+
=
−
−
−
−
x
x
x
x
e
e
y
e
y e
20
2
20
1
20
In the following graph we shows the numerical result obtained with three different integration
methods, Exponential, Runge Kutta of 4
th
order and Euler, using different integration steps
The Expo method reaches the integration length of about x = 19 with a step h = 0.4
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
0
5
10
15
20
Expo method
h = 0.42
The Runge-Kutta method, using a step h = 0.14 gives accurate only for x < 3; after that the
solution begins to diverge.
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
0
5
10
15
20
Runge Kutta 4
h = 0.14
The Euler method, with a step h = 0.11, begins to oscillate even from the first steps. At this step
the Euler method is completely unstable.
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
0
5
10
15
20
Euler method
h = 0.11
As we can see the stability step of the Exponential method is here about 3 times greater than
the others methods.
Documents you may be interested
Documents you may be interested
