Industrial Mathematics - Alexiades
Lab 5
Euler scheme
( Preparation for Project 1 )
Implement the Euler scheme for the scalar IVP:
y′(t) = F(t, y(t)),
t0 ≤ t ≤ tend ,
y(t0) = y0.
The code should read (and print out): the initial point: t0, y0; the final time: tend;
and the number of time-steps to be executed: Nsteps.
The evaluation of F(t,y) should be done in a Function Subprogram FCN( tn, Yn )
which returns the value Fn = F(tn,Yn).
Calculate the time-step Δt = (tend − t0) / Nsteps, implement the Euler time-stepping:
Y0 = y0 ; Yn+1 = Yn + Δt F( tn , Yn ), for n=0,1,...,Nsteps,
with timesteps: tn = t0 + n * Δt , n=0,1,...,Nsteps ,
and output the pairs tn, Yn (and the exact solution when known, see below).
NOTE: No arrays are actually needed...
1. Debug your code on the trivial problem: y'(t) = 2t, 0 ≤ t ≤ 2, y(0) = 1,
by calculating Yexact(tn) = tn*tn + 1, and the error ERRn = ABS(Yexact - Yn).
At each step, output tn, Yn, Yexactn, ERRn,
and calculate the worst overall error: ERRmax = MAX( ERRmax, ERRn).
Output ERRmax at the end of the run.
Before you type in the code, you should write it on paper,
and execute the steps by hand/calculator for, say, Nsteps=4,
to make sure the logic is right. Then enter the code and
check that it finds the same values. Then run it with Nsteps=10.
Then COMMENT OUT PRINTING of tn,Yn, ... and try Nsteps=1000.
Note that since the RHS of this ODE is just a function of t, the solution is simply the integral of F(t)=2t,
so in fact your code is performing numerical integration (usually called "quadrature").
2. Test your code on:
y′(t) = y, 0 ≤ t ≤ 1, y(0) = 1, with Nsteps=10.
Only need to change the FCN, and the exact solution; comment out previous F and enter this one.
Then COMMENT OUT PRINTING of tn,Yn, ... and try Nsteps = 100, 1000, 5000, 10000
to see how the error decreases.
Here the RHS contains the unknown y(t) itself, so this is no longer just a quadrature problem! It is a genuine ODE (separable).
Submit ONLY the following, as "Lab5.txt"
a. ERRmax for Nsteps=10 and for Nsteps=1000 from 1.
=============================================================
b. ERRmax for Nsteps= 1000 and Nsteps=5000 from 2.
=============================================================
c. Your code, including function FCN.
Plotting multicolumn files with gnuplot
You can use gnuplot to plot several curves on the same plot, from columns of a file, as follows:
Gnuplot ignores lines that begin with #, so any lines in
an output file that don't contain plot points
should begin with # (e.g. input values, column labels, etc).
Suppose your Euler code prints output to a file named "out".
It could look something like:
# Input values: t0=0.0000 y0= 1.0000
# tend= 5.000 Nsteps= 10
# tn Yn Yexactn ERRn
0.0 1.0 1.0 0.0
0.2 1.043 1.04 0.3e-02
0.4 1.159 1.16 -0.1e-02
. . .
The pairs (tn, Yn) represent (t, Y(t)) i.e. numerical approximation,
and the pairs (tn, Yexactn) are points of the curve y=Yexact(t).
We can easily plot both curves on the same plot:
gnuplot> plot "out" using 1:2, "out" using 1:3
Better yet, with different colors (lt=linetype, lw=linewidth):
gnuplot> plot "out" u 1:2 lt 1 lw 2, "out" u 1:3 lt 3 lw 2
Linetypes: 1 = red, 2 = green, 3 = blue, 4 = magenta, 7 = black, 0 = black, -1 = solid black
To set labels:
gnuplot> set xlabel "time"
gnuplot> set ylabel "y(t)"
gnuplot> set title "Euler scheme & exact solution: example1"
gnuplot> replot
You can save your settings into a file: gnuplot> save "example1.gplot"
Next time you enter gnuplot, do:
gnuplot> load "example1.gplot"
On linux, you can just do, at the linux prompt: gnuplot < example1.gplot
Better yet, get the template, edit it, and load it.