Matematiska Vetenskaper, Chalmers Tekniska Högskola och Göteborgs Universitet

ALA-C, 2007, studio 2a

Studio 2a

1. In this studio exercise you will use the program mytrig.m to solve some linear systems of ODE's and study their phase portraits. The idea is to compare you numerical computations with the exact solutions that you obtain by using the eigenvalue method. There is also a review on the midpoint rule (cG(1)) (see ALA B) for improved numerical results.
Consider now the eigenvalue problem for the matrix A=[0 1; -1 0].
(i) Find the eigenvalues and corresponding eigenvectors to A. Do this both analytically (by hand) and by using Matlab's [P,D] = eig(A).
(ii) For the same matrix A, solve analytically the initial-value problem
x'(t)= Ax(t), t>0, x(0)=(0,1).
by the eigenvalue method.
(iii) Sketch the phase portrait analytically.
(iv) Now use mytrig.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 the program my_ode2a.m, if you need a hint.

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

4. Try your competence on the following problems:
(a) Consider the following second order ODE:
x''(t) - x(t) = 0, t>0, x(0)=0, x'(0)=1.
(i) Rewrite the problem as a system of two first order equations.
(ii) Solve the problem by hand by using the eigenvalue method.
(iii) Sketch the phase portrait by hand from the information in (ii).
(iv) Use your favourite solver from above (mytrig.m) or the improved solver (my_ode2a.m) to solve the system and plot the phase portrait. By writing a script like

clear
axis([-5 5 -5 5])
hold on
for i=-3:.5:3
for j=-3:.5:3
[t,W]=mytrig([0 2],[i;j],.1);
plot(W(:,1),W(:,2))
pause(0.2)
end
end

you can plot the phase portrait (orbits) for initial values running over the loops in i and j and enjoy how the portrait is slowly coming up.
(b) Repeat the exercise in (a) for the equation
x''(t) + ax'(t) + x(t) = 0, t>0, x(0)=0, x'(0)=1.
for different choices of the constant a.

/nils