Matematik och Datavetenskap, Chalmers Tekniska Högskola och Göteborgs Universitet

ALA-B, 2003, studio 3.2

Home, w 1, w 2, w 3, w 4, w 5, w 6, w 7. Matlab: analysis, linear algebra, facit.

Obligatory work for your individual folder

My ode solver 1

1. Read AMBS Ch 40.4. Write a program based on the forward Euler method for solving the initial value problem for a general system of ordinary differential equations:

u'(t)=f(t,u(t)),     a< t <b,
u(a)=ua.

You must use the template my_ode.m.

Remember that the algorithm is

t0=a
U(t0)=ua

ti=ti-1+h
U(ti)=U(ti-1) + h f(ti-1,U(ti-1))

This program is VERY SIMILAR to my_exp.m. and my_trig.m that you have already written. Use this fact and copy stuff from your old programs.

2. You MUST NOT include any plot commands in the program my_ode. Plotting is a separate task that we may, or may not, want to do after we have solved the differential equation. It is convenient to put such plot commands in a script file called, for example, test1.m:

f='funk';
int=[0 1];
ua=[0;1];
h=1e-3;
[t,U]=my_ode(f,int,ua,h);
figure(1)
plot(t,U) % solution curves
figure(2)
plot(U(:,1),U(:,2)) % phase portrait
axis equal

To execute these commands you type test1 in the matlab command window. This is more convenient than to re-type the commands every time you want to run your program. Note the two kinds of plots: solution curves and phase portrait.

3. Solve the following equations by hand and then use them to test your program. (Note: you have already done this in earlier exercises. Now do it again.)

a.   u'(t)=1/t,    u(1)=0

b.  u'=u,   u(0)=1

c.  u'=-u,   u(0)=1

d.  u''+u=0,   u(0)=0, u'(0)=1

e.  u''-u=0,   u(0)=0, u'(0)=1

You can compare the solution formula and the numerical solution by plotting them in the same figure. For example, if the solution formula is u(t)=3*exp(t):

plot(t,U)
hold on
fplot('3*exp(x)',[0 1],'r')

NOTE: Matlabs built-in function exp.m is NOT AN EXACT SOLUTION, it is a numerical approximation just as the solution that you compute with my_ode and my_exp.

4. The Matlab programs ode23, ode45, ode15s do more or less the same thing as my_ode. Read about it: help ode23 or in helpdesk or in Pärt-Enander-Sjöberg. Use ode23 to solve some of the examples. Compare the numerical solutions obtained by my_ode and ode23.

[t,U]=ode23('funk',int,ua)

5. Documentation. Write down what you have done for your individual folder and show it to your studio teacher before the end of week 4. Include all the test examples carefully written down.

/stig


Last modified: Thu Nov 13 11:54:30 MET 2003