Below is an outline for a possible solution to the project. The computations were done on a sun workstation using MATLAB.

Solution to Assignment One

We run the MATLAB script kepler.m. This script generated the data to answer to two questions.

1. The residuals:

       res_23 = 3.0417e-04
       res_45 = 1.6399e-05
   

   Comparing the norms of the residuals shows that ode45 is ten times more accurate than ode23 for this problem.

2. The number of steps:

       steps_23 = 18
       steps_45 = 61
   

   The length of the vectors on return show that ode45 takes more than 3 times as many steps than ode23 to compute y(2*pi).

Solution to Assignment Two

Using the functions cork.m and cork1.m, the script sol_two.m generates the plots and data to answer the questions of this assignment:

1. Taking t1 = 1.2457, we find x(t1) = 3.7507 with an overshoot error of |x(t1) - L| = 6.8E-4.

2. If we increase L from 3.75 to 3.76, our event happens later, at t2 = 1.2494, with error |x(t2) - L| = 1.6E-3. A increase of L of size 1.0E-2 has led to a increase in t* by |t2 - t1| = 3.7E-3. So the problem is not sensitive to changes in the length L of the cork.

   Looking at the plot of the trajectory, we interpret the slope of x(t) at t=t* as the velocity. As this velocity is larger than one, solving for t in x(t) = L, leads to a minor shift in t as L increases.

3. We take our changes in the initial conditions small enough so we do not have to recalculate the t*.

        x(0) = 0, x'(0) = 1.0e-4: |x(t*) - L| = 9.9861e-04
        x(0) = 1.0e-4, x'(0) = 0: |x(t*) - L| = 5.0083e-04
   

   Looking at the overshoot error, we observe that error is almost twice as much when changing the initial velocity, compared to changing the initial position.

   Comparing the plots of the position x(t) and velocity v(t) around t = 0, we see that v(t) grows faster than x(t). We conclude that the problem is more sensitive to changes in initial velocity than to changes in initial position.

Solution to Assignment Three

We use the script sol_three.m to produce the plots and data.

1. The plots show the solution trajectory to converge to zero, the end vector is endx = [-1.8207e-04 1.3659e-04].

2. When we change the initial conditions from [-4 3] to [-4.001 3.001], the trajectory diverges, ending in endx1 = [5.8134e+04 7.7512e+04]. A tiny change of order 1.0E-3 gives rise to a huge difference in the solution. The solution to this problem is very sensitive to changes in the initial position.

3. The outcome of [v,d] = eig(A) is

   v =
   
     -8.0000e-01   6.0000e-01
      6.0000e-01   8.0000e-01
   
   d =
   
       -1     0
        0     2
   

   Computing c = v\xzero' gives

   c =
   
      5.0000e+00
     -3.5527e-16
   

   which shows that xzero = [-4 3] is five times the first eigenvector with the corresponding negative eigenvalue -1.

   As the general solution to this linear system is a linear combination of the eigenvectors times exp(lambda*t), the solution will diverge when some eigenvalues are positive, unless the starting vector has no component in the direction of eigenvectors with positive eigenvalues.

   Even as the solution to a linear system with positive eigenvalues (or eigenvalues with positive real part) converges, the solution is extremely sensitive to changes in the initial conditions.