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

ALA-B, 2003, studio 4.2

My ode solver 2

1. Read AMBS Ch 40.7. In this studio exercise you will modify your program my_ode.m so that it uses the midpoint method instead of the Euler forward method.

First use my_ode.m to solve the system

w'=Aw,    a<t<b,
w(0)=w0

where w=[u; u'], w0=[0; 1], A=[0 1; -1 0]. Use a rather large step (h=1e-1).

Plot the solution

figure(1)
plot(t,W) % solution curves

figure(2)
plot(W(:,1),W(:,2)) % phase plane plot
axis equal

Note that the amplitude is increasing. The reason for this is that the Euler forward method does not conserve the energy in this equation. Instead the energy is increasing.

2. The Euler forward method solves the initial value problem for a general system of ordinary differential equations:

u'(t)=f(t,u(t)),     a< t <b,
u(a)=u_a,

by the algorithm

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

The midpoint method is

U(t_i) =U(t_i-1)+h f( ( t_i-1+t_i)/2, ( U(t_i-1)+U(t_i) )/2)

Note that this is a fixed point equation for V=U(t_i), of the form V=g(V):

V =U(t_i-1)+h f( ( t_i-1+t_i)/2, ( U(t_i-1)+V )/2)

This equation has a unique solution if the Lipschitz constant of g is smaller than 1, L_g = h L_f/2 < 1, that is, if h < 2/ L_f.

Implement the corresponding fixed point iteration in a new program my_ode2.m:

V_k =U(t_i-1)+h f( ( t_i-1+t_i)/2, ( U(t_i-1)+V_k-1 )/2)

Use the tolerance for the iteration tol=h^3. See my program my_ode2.m, if you need a hint.

3. Now solve the system for the trigonometric functions again. Does the amplitude increase?

4. Advanced.
a. Write a program for the Euler backward method. See my program my_ode1.m, if you need a hint.
b. Use Newton's method instead of fixed point iteration.
c. Choose the step adaptively.
Hints: see the ode programs in facit.

/stig