Industrial Mathematics - Alexiades
Lab 6
Euler scheme for systems of ODEs
(preparation for project 2)
It is easy to modify your Euler scheme code for a system of (two) ODEs:
y1′(t) = F1(t, y1(t), y2(t)),
y1(t0) = y10
y2′(t) = F2(t, y1(t), y2(t)),
y2(t0) = y20
Copy the file with your scalar Euler scheme code to a file with a new name
(e.g. Euler2.m), and make the necessary changes.
You can simply rename the single old unknown to Y1n, and introduce
the 2nd unknown Y2n (and y20, and F2),
and a 2nd Euler update for Y2.
The functions F1 and F2 should be evaluated in one FCN subprogram.
Debug your code on the easy problem, for 0 ≤ t ≤ 1:
y1′ = y2 , y1(0) = 1,
y2′ = y1 , y2(0) = 0.
This can be solved exactly: Setting z=y1, we have z'= y1' = y2, and z''= y2'= y1 = z,
so this system is equivalent to the 2nd order linear ODE
z'' − z = 0, with ICs: z(0)=1, z'(0)=0.
Solve it to show that the unique solution is:
y1(t) ≡ z(t) = cosh(t) ,
y2(t) ≡ z'(t) = sinh(t)
For debugging, calculate Y1exact(tn) = cosh(tn)
and the error ERRn = ABS(Y1exact - Y1n).
At each time step, output:
tn Y1n Y1exactn ERRn
and calculate the worst overall error:
ERRmax = MAX( ERRmax, ERRn).
Output ERRmax at the end.
After debugging, comment out the printing!
Make runs with N=1000 , N=5000 , N=10000 time-steps.
Observe how the error decreases.
Submit ONLY the following, in a text file Lab6.txt:
-----------------------------------------------------------
Name, date, Lab6
-----------------------------------------------------------
from run with N=1000: input , y1(1) , ERRmax
-----------------------------------------------------------
from run with N=5000: input , y1(1) , ERRmax
-----------------------------------------------------------
your code, and your FCN.