MCS 471 Project Five Solution

Below is an outline for the solution of project five.

1. the scripts

  1. assignment_one.m 
    % solution to assignment one:
%
% 1. first solve the problem
%
kepler = inline('1/(1-.2*cos(y))','x','y');
xspan = [0 2*pi];
yzero = 0;
[x y] = ode45(kepler,xspan,yzero);
z = y - .2*sin(y) - x;
fprintf('largest value from the equation : %e\n',max(abs(z)));
%
% 2. examine the sensitivity to initial value
%
[x0 y0] = ode45(kepler,xspan,yzero+1.0);
[x1 y1] = ode45(kepler,xspan,yzero+1.0e-1);
[x2 y2] = ode45(kepler,xspan,yzero+1.0e-2);
[x3 y3] = ode45(kepler,xspan,yzero+1.0e-3);
figure
plot(x,y)
hold on
plot(x0,y0)
plot(x1,y1)
plot(x2,y2)
plot(x3,y3)
fprintf('error from 1.0e-1 change: %e\n',abs(y0(end)-y(end)));
fprintf('error from 1.0e-1 change: %e\n',abs(y1(end)-y(end)));
fprintf('error from 1.0e-2 change: %e\n',abs(y2(end)-y(end)));
fprintf('error from 1.0e-3 change: %e\n',abs(y3(end)-y(end)));
  2. assignment_two.m 
    % solution to assignment two:
%
%  1. first solve with damping for k = 0.2
%
f = inline('[y(2); -0.2*y(2)-sin(y(1))]','t','y');
yzero = [1 0];                          
tspan = [0 50]; 
[t y] = ode45(f,tspan,yzero); 
figure
plot(t,y(:,1)) 
fprintf('value for theta at the end: %e\n',abs(y(end)));
[t1 y1] = ode45(f,tspan,yzero+0.01);
err = abs(y1(end,1)-y(end,1));
fprintf('error at end with 0.01 change in theta(0): %e\n',err);
%
%  2. solve this problem with k = -0.04
%
f = inline('[y(2); 0.04*y(2)-sin(y(1))]','t','y');
yzero = [1 0];                          
tspan = [0 60]; 
[t,y] = ode45(f,tspan,yzero); 
figure
plot(t,y(:,1)) 
fprintf('value for theta at the end: %e\n',abs(y(end)));
[t1 y1] = ode45(f,tspan,yzero+0.01);
err = abs(y1(end,1)-y(end,1));
fprintf('error at end with 0.01 change in theta(0): %e\n',err);
  3. assignment_three.m 
    % solution to assignment three:
%
%  1. first solve with ode45
%
f = inline('[1195*y(1) - 1995*y(2); 1197*y(1) - 1997*y(2)]','t','y');
yzero = [2 -2];                          
tspan = [0 1]; 
[t y] = ode45(f,tspan,yzero); 
subplot(2,1,1);
plot(t,y(:,1)) 
hold on
x1 = 10*exp(-2*t) - 8*exp(-800*t);
plot(t,x1);
subplot(2,1,2);
plot(t,y(:,2)) 
x2 = 6*exp(-2*t) - 8*exp(-800*t);
hold on
plot(t,x2);
fprintf('maximal error for x : %e\n',max(abs(x1-y(:,1))));
fprintf('maximal error for y : %e\n',max(abs(x2-y(:,2))));
fprintf('number of steps : %d\n',length(t));
%
%  2. now solve with ode15s
%
[t y] = ode15s(f,tspan,yzero); 
subplot(2,1,1);
plot(t,y(:,1)) 
hold on
x1 = 10*exp(-2*t) - 8*exp(-800*t);
plot(t,x1);
subplot(2,1,2);
plot(t,y(:,2)) 
x2 = 6*exp(-2*t) - 8*exp(-800*t);
hold on
plot(t,x2);
fprintf('maximal error for x : %e\n',max(abs(x1-y(:,1))));
fprintf('maximal error for y : %e\n',max(abs(x2-y(:,2))));
fprintf('number of steps : %d\n',length(t));

2. the main MATLAB session

We executed the three scripts in a MATLAB session, after typing diary proj5sol, the content of the file proj5sol is 
assignment_one
largest value from the equation : 7.215717e-06
error from 1.0e-1 change: 9.999998e-01
error from 1.0e-1 change: 1.000000e-01
error from 1.0e-2 change: 1.000043e-02
error from 1.0e-3 change: 9.999986e-04
assignment_two
value for theta at the end: 4.881709e-03
error at end with 0.01 change in theta(0): 5.830269e-05
value for theta at the end: 2.415193e+00
error at end with 0.01 change in theta(0): 1.164883e+00
assignment_three
maximal error for x : 3.839707e-03
maximal error for y : 3.839707e-03
number of steps : 997
maximal error for x : 2.349006e-03
maximal error for y : 2.349014e-03
number of steps : 59
diary off

3. the plots

  1. first plot
  2. second plot
  3. third plot

4. answers to the questions

  1. Assignment One

    From plugging the (x,y) tuple into the equation, we see that the residual is about 10^(-6). An error on the initial value of the order 10^(-k) results in an error at the end of the solution trajectory of the same magnitude, see the table: 

          error from 1.0e-1 change: 9.999998e-01
      error from 1.0e-1 change: 1.000000e-01
      error from 1.0e-2 change: 1.000043e-02
      error from 1.0e-3 change: 9.999986e-04
    So the problem is well conditioned. From the first plot we see the parallel trajectories.

  2. Assignment Two

    The second plot above shows the solution trajectory, adjusted with a time span to 50 to see the convergence to an equilibrium. In the third plot we see the divergence, when the pendulum swings out of control for k = -0.04.

    The sensitivity analysis shows that the converging pendulum is well conditioned (error of 5.8e-5), while the diverging pendulum is very sensitive to changes in the initial values (error of 1.1e+0). We can interpret this from the plots: decreasing amplitudes will decrease errors and increasing amplitudes will increase errors.

  3. Assignment Three

    Both ode45 and ode15s give the same accuracy, but ode45 needs 997 steps, while ode15s needs only 57 steps to reach a precision of 1.0E-3.

5. downloading the scripts