MATLAB Assignment 1

Math 222 Spring 2018

You learned in class the Euler method for approximating a solution to the differential equation

$$ \frac{d y}{d t} = f(t,y). $$

Given an initial condition $y(0)=y_0$, the method gives a sequence of approximations $y_n \approx y(t_n)$ using the recursion scheme

$$ y_{n+1} = y_n + h f(t_n, y_n). $$

where $t_n =  n h$ and $h$ is the time step chosen by the user.

Contents

Warm up

A simpler recursion relation would look like this

$$ y_n = \sum_{j=1}^{n} j$$

N = 10;
y = zeros(1,N+1);
y(1) = 0;
for j=1:N
    y(j+1) = y(j) + j;
end
plot(0:N,y,'o')
xlabel('n');ylabel('y_n')

The Euler method

Let's solve the problem

$$ \frac{dy}{dt} = \sin{t} y; y(0) = 1.$$

This has exact solution

$$ y = e^{1-\cos{t}}.$$

tFinal = 2*pi;
N = 100;
h=tFinal/N;
t=linspace(0,tFinal,N+1); % type 'help linspace' to see how this works
y=zeros(1,N+1);
yExact=exp(1-cos(t));
y(1) = 1; % setting the initial condition
for n=1:N
    y(n+1) = y(n) + h * sin(t(n))*y(n);
end
plot(t,y,t,yExact,'--'); xlabel('t'); ylabel('y'); title('Look, ma! I solved an ODE in MATLAB!');

error100= abs(y(N+1)-yExact(N+1));

That doesn't do so well, let's double the number of points and try again

N = 200;
h=tFinal/N;
t=linspace(0,tFinal,N+1); % type 'help linspace' to see how this works
y=zeros(1,N+1);
yExact=exp(1-cos(t));
y(1) = 1; % setting the initial condition
for n=1:N
    y(n+1) = y(n) + h * sin(t(n))*y(n);
end
plot(t,y,t,yExact,'--')
xlabel('t'); ylabel('y'); title('Look, ma! I solved another ODE in MATLAB!');

error200= abs(y(N+1)-yExact(N+1));

fprintf('The error went down by a factor of  %f.\n',error100/error200);
The error went down by a factor of  1.953899.

Your assignment

Question 1: Apply the Euler method to solve

$$ \frac{dy}{dt} =y(3 - t y); y(0) = 1$$

from $t=0$ to $t= 2$.

Cut and paste the above snippet of code into the MATLAB editor. Modify as necessary and save it as a .m file in order to do the following:

Note that this problem has exact solution

$$ y(t) = \frac{9}{3 t- 1 +10 e^{-3 t}}. $$

You should hand in a printout of your code, any figures you print, and a report of a few sentences answering all the questions.

Bonus Question 2 for the brave The Euler method is an example of a first order method. This means that the error is proportional to $h^2$, so that doubling the number of points roughly cuts the error roughly by a factor of four.

The improved Euler method is what's known as a two-step method. It requires two evaluations of the function per time step. At each step, we define two extra quantities before computing $y_{n+1}$

$$ k_1= f(t_n,y_n) $$

$$ k_2 = f(t_{n+1}, y_n + h k_1) $$

$$ y_{n+1} = y_n + \frac{h}{2}(k_1 + k_2).$$

Program this method to solve the same problem from Question 1 and