Math 344 - Numerical Methods - Homework Assignment 9

Homework Assignment 9 -- Turn in on March 17

Review the Euler Method, the modified Euler Method, Heun Method, and the Midpoint Method given in the last homework assignment. The last three methods are also called the 2nd order Runge-Kutta methods. Review also the 4th order Runge-Kutta methods and particularly the one called the standard 4th order Runge-Kutta method.

Here are function programs in MatLab for these methods.

Euler Method:

    % eulfun.m
    % Euler's Method
    %
    % yv=eulfun(xv,y0,h,n);
    %
    function yv=eulfun(xv,y0,h,n);
    yv(1)=y0;
    for ii=2:n+1;
        f1=fun(xv(ii-1),yv(ii-1));
        yv(ii)=yv(ii-1)+h*f1;
    end;

Heun's Method:

    % heun.m
    % heun -- Heun's Method
    %
    % yv=heun(xv,y0,h,n);
    %
    function yv=heun(xv,y0,h,n);
    yv(1)=y0;
    for ii=2:n+1;
        f1=fun(xv(ii-1),yv(ii-1));
        tempy=yv(ii-1)+2*h/3*f1;
        f2=fun(xv(ii-1)+2*h/3,tempy);
        yv(ii)=yv(ii-1)+h*(f1+3*f2)/4;
    end;

Modified Euler Method:

    % modeul.m
    % Modified Euler's Method
    %
    % yv=modeul(xv,y0,h,n);
    %
    function yv=modeul(xv,y0,h,n);
    yv(1)=y0;
    for ii=2:n+1;
        f1=fun(xv(ii-1),yv(ii-1));
        tempy=yv(ii-1)+h*f1;
        f2=fun(xv(ii-1)+h,tempy);
        yv(ii)=yv(ii-1)+h*.5*(f1+f2);
    end;

Midpoint Method:

    % midpt.m
    % Midpoint Method
    %
    % yv=midpt(xv,y0,h,n);
    %
    function yv=midpt(xv,y0,h,n);
    yv(1)=y0;
    for ii=2:n+1;
        f1=fun(xv(ii-1),yv(ii-1));
        tempy=yv(ii-1)+h/2*f1;
        f2=fun(xv(ii-1)+h/2,tempy);
        yv(ii)=yv(ii-1)+h*f2;
    end;

the Standard Runge-Kutta Method:

    % rk4.m
    % rk4 -- a 4th order Runge Kutta Method
    %
    % yv=rk4(xv,y0,h,n)
    %
    function yv=rk4(xv,y0,h,n);
    yv(1)=y0;
    for ii=2:n+1;
        f1=fun(xv(ii-1),yv(ii-1));
        tempy=yv(ii-1)+h/2*f1;
        f2=fun(xv(ii-1)+h/2,tempy);
        tempy=yv(ii-1)+h/2*f2;
        f3=fun(xv(ii-1)+h/2,tempy);
        tempy=yv(ii-1)+h*f3;
        f4=fun(xv(ii-1)+h,tempy);
        yv(ii)=yv(ii-1)+h*(f1+2*f2+2*f3+f4)/6;
    end;

In the following example, we will use these 5 methods to approximate the solution of the initial value problem

    y'=t+y,  y(1)=2, for t in [1,5].

Since we know the true solution is: y(t)=-(t+1)+4e^t-1, the true error for each method at each t_k is also computed. Both approximation y₁, y₂, y₃, ... and error e₁, e₂, e₃, ... where e_k = |y(t_k)-y_k| are plotted. Here are the MatLab commands used for this example and plots (h=0.2) showed on the screen.

    %
    % Solve y'=t+y, y(1)=2, t in [1,5]
    %
    clear
    a=1;
    b=5;
    y0=2;
    h=input('stepsize h = ');
    tv=a:h:b;
    ysol=-(tv+1)+4*exp(tv-1);
    n=length(tv);
    yeuler=eulfun(tv,y0,h,n);
    yheun=heun(tv,y0,h,n);
    ymodeul=modeul(tv,y0,h,n);
    ymidpt=midpt(tv,y0,h,n);
    yrk4=rk4(tv,y0,h,n);
    clf
    subplot(211),plot(tv(1:n),yeuler(1:n),':');
    hold
    plot(tv(1:n),yheun(1:n),'-.');
    plot(tv(1:n),ymodeul(1:n),'--');
    plot(tv(1:n),ymidpt(1:n),'-');
    plot(tv(1:n),yrk4(1:n),':');
    plot(tv(1:n),yrk4(1:n),'o');
    plot(tv(1:n),ysol(1:n),':');
    plot(tv(1:n),ysol(1:n),'*');
    title('One Step Methods: dy/dt=t+y, y(1)=2, t in [1,5]')
    text(1.5,200,'... Euler Method')
    text(1.5,180,'-.- the 2nd order Runge-Kutta Methods')
    text(1.5,160,'.o. the Std 4th order Runge-Kutta Method')
    text(1.5,140,'.x. the solution y(t)')
    hold off
    erreul=abs(ysol(1:n)-yeuler(1:n));
    errheun=abs(ysol(1:n)-yheun(1:n));
    errmode=abs(ysol(1:n)-ymodeul(1:n));
    errmidpt=abs(ysol(1:n)-ymidpt(1:n));
    errrk4=abs(ysol(1:n)-yrk4(1:n));
    subplot(212),plot(tv(1:n),erreul(1:n),':')
    hold
    plot(tv(1:n),errheun(1:n),'-.')
    plot(tv(1:n),errmode(1:n),'--')
    plot(tv(1:n),errmidpt(1:n),'-')
    plot(tv(1:n),errrk4(1:n),':')
    plot(tv(1:n),errrk4(1:n),'o')
    title('Errors for each method: |y(t_k)-y_k|')
    text(1.5,60,'... Euler Method')
    text(1.5,50,'-.- the 2nd order Runge-Kutta Methods')
    text(1.5,40,'.o. the Std 4th order Runge-Kutta Method')
    hold off

The graphs are given at the end of this assignment.

Use the Euler Method, the modified Euler Method, Heun Method, the Midpoint Method, and the standard 4th order Runge-Kutta method to approximate numerically the solutions of the following two initial value problems with given h. Plot the sequences of y_k and the true errors if the true solution is known.
- Consider the initial value problem:
```
    y'=ty+2t-t³, y(0)=0,  for x in [0,2]. 
```
  Use h=0.2, h=0.1 and h=0.01. Know that y(t)=t².
- Consider the initial value problem:
```
    y'=t²+y², y(0)=1, for x in [0,0.9].
```
  Use h=0.1, h=0.05 and h=0.01. The true solution is not known. Use the approximation generated by the standard 4th order Runge-Kutta method as the solution to compute the errors of other methods.
Page 462, #6, #7.
Give the formula for y_k+1 when the 3rd Taylor polynomial is used to approximate the solution y(t). Then compute y₁ for the initial value problem:
```
    y'=2ty², y(0)=1, h=0.2.
```
The true solution is: y(t)=1/(1-t²). Find the true error for y₁.

Here are the graphs: