ESAM 446-1

Problem Set 1

Due October 2 in class.

1. ODE: Nonlinear Duffing Oscillator

Consider the ordinary differential equation describing the frictionless motion of a particle in a quartic potential,

d2 u dt2 = -¶u V(u), with V(u) = - 1 2 u2+ 1 4 u4.

(1)

Write a MATLAB program to solve this equation numerically. To monitor the quality of your simulation have the program also calculate the total energy E = [1/2] (¶_tu)² + V(u) of the system. For all cases to be studied use the initial condition u(t = 0) = 1.414214, ¶_t u(t = 0) = 0, which explores a sensitive aspect of the system, because the energy is just large enough for the particle to cross the maximum of the potential at u = 0.

For cases 1. and 3. (below) first validate your code by establishing that it gives the correct order of convergence (order of the scheme) by calculating the solution for t = 4. Determine the dependence of the change in total energy E(t = 4)-E(t = 0) as a function of the time step and plot it double-logarithmically. Also plot double-logarithmically the difference between the values for u(t = 4) in successive approximations as a function of the time step.

After that, in all 4 cases below run the code all the way to t = 100 and study that result.

1. Familiarize yourself with the provided template and use it for some preliminary runs. It uses only simple forward Euler. Can you choose the time step Dt small enough to keep the (numerically induced) energy change below |DE| = 0.1? How does the change in energy manifest itself in the solution?

2. Improve the program by making the time step adaptive and by using Richardson extrapolation to control the time step as well as increase the order of the scheme. Modify stepeuler.m for these modifications. Control the time step such that the difference between the approximations at Dt and Dt/2 varies at most by a factor of 2 over the whole run. Can you obtain qualitatively the correct behavior for 0 < t < 100? How many time steps are needed to get the (numerically induced) change in total energy below |DE| = 10^-7? Does the final value u(t = 100) show already the correct scaling with Dt in this regime?

3. Replace the fixed-time-step Euler code by a fixed-time-step fourth-order Runge-Kutta code. How many steps do you need to reproduce respectively the accuracy achieved with the two versions of your Euler code? Can you get u(t = 100) to be accurate within 1%?

4. Replace the variable-time-step Euler scheme by a variable-step fourth-order Runge-Kutta scheme. How do you have to change the criteria for the modification of Dt as compared to the corresponding Euler code? How many steps do you need to get u(t = 100) accurate to within 1%? To compare the Euler and the Runge-Kutta codes it is better to compare the number of floating point operations rather than time steps (use the command flops for that).

• Write the second-order equation as two coupled first-order equations and use the forward Euler and Runge-Kutta scheme as given in class.

• When decreasing the time step always decrease it by a factor of 2.

• For the log-log plots use the command loglog instead of plot. Do not take the logarithm of the data and plot that result using plot.

• You don't have to hand in your code. Your write-up should demonstrate that your code is working properly, for instance by demonstrating the correct scaling of the error of your code with the time step (see above).

• Hand in:

  • appropriately labeled plots of the convergence tests (energy difference and u at the respective final time as a function of time step)
  • answers to the explicitly stated questions supported by plots and/or tables.
  • appropriately labeled plots of your other measurements/results as needed to support your statements and conclusions
  • discussion of your results with comparison between the schemes and any other observations you made