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.

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").
What quadrature rule does it amount to ??? (rectange? trapezoidal? Simpson?)

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 function y(t) itself, so this is no longer just a quadrature problem!
This example is discussed in the CSEP e-book, chapter: ODE, 2.2.1 Taylor series methods, Euler method (Example 8).

Submit ONLY the following:
  a. ERRmax for Nsteps=10 and for Nsteps=1000   from 1.
  =============================================================
  b. ERRmax for Nsteps= 1000 and Nsteps=5000   from 2.
  =============================================================
  c. Your code.



            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    1.0      1.0         .0
       .2    1.043    1.04        .3e-02
       .4    1.159    1.16       -.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" using 1:2 lt 1 lw 2, "out" using 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 in 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.