# Homework 6, due 08.00, 29/10-2014¶

## Numerical solution of ODEs¶

In this homework you will try out some numerical methods for solving ODEs.

## Duffing’s attractor¶

This problem is intended to demonstrate that there is a difference between good and not so good numerical methods. Consider the second order ODE:

Your assignment is to compute the solution at \(t=70\) as accurately as your computer allows. First try forward Euler:

then the fourth order accurate Runge Kutta method:

Here \(\dot{x} = f(t,x(t))\) is a system of first order ODEs so you should first rewrite the second order ODE that form by introducing \(x_1=z, \, x_2=\dot{z}\).

Use the values

as “exact” and plot the error \(e=\sqrt{(x_1-x^{ex}_1)^2+(x_2-x^{ex}_2)^2}\) as a function of \(\Delta t\). If \(\Delta t\) is small enough (inside the asymptotic regime) the error is expected to follow a power law \(e \approx C (\Delta t)^p\). Determine \(C\) and \(p\) from your plot and estimate the number of time steps required to reach \(e < 10^{-5}\) with forward Euler. Is your computer fast enough?

**Note:** Unless you take \(\Delta t\) small enough the above
methods will not be stable. Take \(\Delta t < 8.75\cdot10^{-3}\)
for Euler and \(\Delta t < 0.35\) for Runge-Kutta to be on the
safe side.

Report the order of convergence by plotting the error vs. the timestep
for both method. Also plot the trajectory of the solution,
`plot(x(:,1),x(:,2))` for some different timesteps.

## The Wave equation in one dimension¶

Consider the wave equation in one dimension with a variable coefficient

with homogenous boundary conditions

In this homework we discretize the first order system

as described in the file `/codes/fortran/HWK6/rhside.f90`. Note that
as both \(u\) and \(v\) are zero on the boundary we only have
to evolve the solution on the interior gridpoints \(x_i, \, i =
1,\ldots,n_x-1\). After discretization in space we thus have to solve a
system of ODE of size \(2(n_x-1)\)

Here we have used the notation

See also `rhside.f90`.

Start off by compiling and running the program
`/codes/fortran/HWK6/homework6.f90`. The program currently uses
second order Runge-Kutta, or Heun’s method to solve the above system
of ODE. This method is not stable but it takes a little while before
the solution bows up. Run the code with `nx = 100` and `nx = 200`
and look at the output by the Matlab program `plotresults.m`. What
if you increase `nx` further? When does the solution blow up? after
a fixed number of timesteps, at a fixed time? Keep refining and try to
find a pattern (for example by plotting the timestep the solution
reaches one million as a function of number of gridpoints).

- Modify the code to use Runge-Kutta 4 and try some different values of
amp. What happens if \(a(x) < 0\)?- Try and change \(a(x)\), can you make the waves get “trapped”?
- For the case with
amp = 0.d0, slowly increasecfluntil the solution blows up. How can we relate the time-stepping stability of this problem to the stability of the Dahlquist equation \(y' = \lambda y\)? As we saw in class, for a linear system of equations \({\bf y}' = {\bf A}{\bf y}\) the eigenvalues of the matrix plays the role of \(\lambda\) in the scalar case. But what is the matrix in this case? And what are the eigenvalues? The right hand side is in fact nothing but a linear transformation so (consult your linear algebra book if you for got this) we may reconstruct the matrix, column by column if we plug in unit vectors info the linear transformation. Look at the program/codes/fortran/HWK6/getmatrix.f90to see how we can construct the matrix. The program also computes the eigenvalues. Look up the documentation of the Lapack subroutinedgeevto see what is being output in the fileslam_re.txtandlam_im.txt.- Find the stability domain of Runge-Kutta 4 in a numerical analysis book or online. Compute the eigenvalues for a couple of different values of
hx(chosen by settingnx). Multiply the eigenvalues bydt = cfl*hx. How large cancflbe if all of the \(\Delta t \lambda\) are to remain iside the stability domiain of Runge-Kutta 4? Does this agree with your computational / experimental results from 3.?- Finally, as you now know the matrix, implement Backward Euler as an option for time-stepping in your program. Let \({\bf W} = [u\, v]^T\) Then you can write backward Euler \({\bf (I }-\Delta t {\bf A) W(t_{n+1})}= {\bf W(t_{n})}\). Read the documentation for
DGETRFandDGETRSand use them to factor (once before the time-stepping loop) and solve to get the next timestep. Is is true that you can take any size timesteps? Compute the solution to some problems using RK4 and backward Euler and compare the accuracy and computational time between them.