AMSC 661 - Homework 3

Daniel Dunlavy  Darran Furnival  John Harlim

April 18, 2002

Throughout the project we use the following Fourier transforms:

u_{j}^{n} = u(x_{j}, t_{n}) = ∑_{j=1}^{N} û_{k}^{n}e^{ikx_{j}}  j = 1, 2, …, N

û_{k}^{n} = û_{k}(t_{n}) = \frac{1}{N} ∑_{j=1}^{N} u_{j}^{n}e^{−ikx_{j}}  k = 1, 2, …, N

Transforms for f_{j}^{n} and f̂_{k}^{n} are defined in a similar way.

(1) Explore spectral method and nonlinear wave equations.

In question 1 we use three numerical methods to solve the ODE’s that result from taking Fourier transforms of the given PDE’s, namely the leap frog method, the second order leap frog method and a fourth order Runge-Kutta method. We here look at the absolute stability of these methods.

To examine absolute stability we apply each method to the model problem

u_{t} = λu

and then look at the resulting stability polynomial. A method will be absolutely stable if and only if the roots of the stability polynomial lie in the unit disk in the complex plane with those roots on the boundary of the disk being simple.

Applying the leap frog formula to the model problem gives

\frac{u^{n+1} − u^{n−1}}{2Δt} = λu^{n}

⇔ u^{n+1} − 2Δtλu^{n} − u^{n−1} = 0.

Setting σ = Δtλ then gives

(u^{2} − 2σu − 1)u^{n} = 0

where we have shifted the index. The stability polynomial is then given by

π(u) = u^{2} − 2σu − 1.

We can now use the following theorem to examine the roots of the stability polynomial.

Theorem 1 The roots x_{1,2} of ax^{2} + 2bx + c = 0 satisfy |x_{1,2}| ≤ 1 if and only if

• |c| < |a| and 2|āb − b̄c| ≤ |a|^{2} − |c|^{2}

or

• |c| = |a|, āb = b̄c and |b| ≤ |a|.

Applying these conditions to the stability polynomial of the leap frog method we see that the first set of conditions are violated but that the second set are fulfilled provided that σ is pure imaginary and that |σ| ≤ 1. Therefore the leap frog method will be stable if σ ∈ (−i, i); we neglect the end points of the interval as they will make little practical difference and it saves us worrying about such roots not being simple.

For the second order leap frog method the model problem becomes

u_{tt} = λu

which becomes

\frac{u^{n+1} − 2u^{n} + u^{n−1}}{Δt^{2}} = λu^{n}

⇔ u^{n+1} − (2 + λΔt^{2})u^{n} + u^{n−1} = 0.

Setting σ = λΔt^{2} gives, with a shift of index

(u^{2} − (2 + σ)u + 1)u^{n} = 0

from which we see that the stability polynomial is

π(u) = u^{2} − (2 + σ)u + 1.

Applying the above theorem we readily ascertain that σ must be real and \frac{1}{2}|2 + σ| ≤ 1 for the roots to lie in the unit disk on the complex plane. Therefore we will get absolute stability provided that σ ∈ (−4, 0) where we drop the end points for the same reason as above.

The Runge-Kutta forth order method when applied to u_{t} = f(u, t) is given by

a = Δtf(u^{n}, t_{n}),

b = Δtf(u^{n} + a/2, t_{n} + Δt/2), c = Δtf(u^{n} + b/2, t_{n} + Δt/2), d = Δtf(u^{n} + c, t_{n} + Δt),

u^{n+1} = u^{n} + \frac{1}{6}(a + 2b + 2c + d).

Applying this to the model problem and setting σ = λΔt gives

u^{n+1} = u^{n} + \frac{1}{6}(σu^{n} + 2σ(u^{n} + \frac{σ}{2} u^{n}) + 2σ(u^{n} + \frac{σ}{2} (u^{n} + \frac{σ}{2} u^{n}))+

σ(u^{n} + σ(u^{n} + \frac{σ}{2} (u^{n} + \frac{σ}{2} u^{n}))))

Figure 1: Stability region for forth order Runge-Kutta method.

The stability polynomial can be arrived at by dividing this through by u^{n} so we are interested in looking at the roots of

u = 1 + \frac{1}{6}(σ1 + 2σ(1 + \frac{σ}{2} ) + 2σ(1 + \frac{σ}{2} (1 + \frac{σ}{2} 1))+

σ(1 + σ(1 + \frac{σ}{2} (1 + \frac{σ}{2} ))))

which simplifies to

u = 1 + σ + \frac{1}{2}σ^{2} + \frac{1}{6}σ^{3} + \frac{1}{24}σ^{4}.

The roots of this will satisfy

1 + σ + \frac{1}{2}σ^{2} + \frac{1}{6}σ^{3} + \frac{1}{24}σ^{4} = 0

and we then want for stability

|1 + σ + \frac{1}{2}σ^{2} + \frac{1}{6}σ^{3} + \frac{1}{24}σ^{4}| ≤ 1.

The boundary of this region is shown plotted in Figure 1. An argument from complex analysis can then be used to show that the region of absolute stability is the region lying inside this curve.

(a) 1st order variable coefficient wave equation

u_{t} + c(x)u_{x} = 0 \{ \begin{matrix}c(x) = .2 + sin^{2}(x + 2) \\ u(x, 0) = exp(−100(x + 2)^{2}) \end{matrix}

We assume periodic boundary conditions for x ∈ [−π, π] and use the following spatial grid:

x_{j} = −π + \frac{2π}{N} j, j = 1, …, N .

We then apply the Fourier transform to the discretized PDE to get a system of ODE’s:

\frac{d}{dt}û_{k} + f̂_{k} = 0 ,

where f(x, t) = c(x)u_{x}(x, t). In order to compute f̂_{k}, we note that partial derivatives of u with respect to x are easy to compute in the Fourier domain:

\{(ik)^{m}û_{k}\} ⟶^{fft} \{(\frac{∂^{m}}{∂x^{m}} u)_{j}\}  (1)

Using the leap frog method, we thus have the following iteration

\frac{û_{k}^{n+1} − û_{k}^{n−1}}{2Δt} + f̂_{k}^{n} = 0 ,

with û_{k}^{0} = F(u(x, 0)), the Fourier transform of the given initial condition. To start this implicit method, we assume a constant wave speed of c = .2 and extrapolate back in time. For the purposes of stability analysis we look at

u_{t} + c_{max}u_{x} = 0
⇔ u_{t} + \frac{6}{5}u_{x} = 0.

In the Fourier domain this becomes

û_{t} + ik \frac{6}{5}û = 0.

From the general stability analysis for the leap frog method we require λΔt ∈ (−i, i). We can see that here λ = −6ik/5 so the stability condition becomes

|λΔt| < 1 ⇒ Δt < \frac{5}{6} \frac{1}{|k|} ≤ \frac{5}{3N}

where we have taken the maximum value of k to be N/2. This is seen to be in good agreement with the empirical results where the method became unstable at Δt = 1.83N^{−1}. Figure 1(a) presents plots of the solution evolving over time for values of Δt giving stable and unstable results.

\begin{matrix}Δt = 1.82/N \\ Δt = 1.83/N \end{matrix}

Figure 1(a): u_{t} + c(x)u_{x} = 0 with N = 128 (b) 2nd order wave equation

u_{tt} = u_{xx} \{ \begin{matrix}u(x, 0) = exp(−10(x)^{2}) \\ u_{t}(x, 0) = 0 \end{matrix}

We assume periodic boundary conditions for x ∈ [−π, π] and use the same spatial grid as in 1(a). Applying the Fourier transform to the discretized PDE, we get the following system of ODE’s:

\frac{d^{2}}{dt^{2}} û_{k} = −k^{2}û_{k} ,

where we computed the partial derivatives in the Fourier domain, using (1).

Using the second order leap frog method, we have the following iteration

\frac{û_{k}^{n+1} − 2û_{k}^{n} + û_{k}^{n−1}}{(Δt)^{2}} = −k^{2}û_{k}^{n}

with û_{k}^{0} = F(u(x, 0)), the Fourier transform of the given initial condition. To start this implicit method, we use the initial condition u_{t}(x, 0) = 0 and set û_{k}^{−1} = û_{k}^{0}, since u is initially not changing with respect to t.

For the second order leap frog we want λΔt^{2} ∈ (−4, 0). We see that in the Fourier domain we have λ = −k^{2} so for stability we require

Δt^{2} < \frac{4}{k^{2}} ≤ \frac{16}{N^{2}} ⇒ Δt < \frac{4}{N} .

This is in perfect agreement with the results obtained empirically.

Figure 1(b) presents plots of the solution evolving over time for values of Δt giving stable and unstable results.

\begin{matrix}Δt = 3.99/N \\ Δt = 4.00/N \end{matrix}

Figure 1(b): u_{tt} = u_{xx} with N = 80 (c) KdV equation and solitary waves

u_{t}+uu_{x}+u_{xxx} = 0  \{ \begin{matrix}u(x, 0) = 3A^{2} sech^{2}(A(x + 2)/2) + 3B^{2} sech^{2}(B(x + 1)/2) \\ A = 25 \\ B = 16 \end{matrix}

We assume periodic boundary conditions for x ∈ [−π, π] and use the same spatial grid as in 1(a). Applying the Fourier transform to the discretized PDE, we get the following system of ODE’s:

\frac{d}{dt}û_{k} + f̂_{k} − ik^{3}û_{k} = 0 ,

where f(x, t) = u(x, t)u_{x}(x, t) and the partial derivatives are computed in the Fourier domain, using (1). Since this is a stiff ODE, we use an integrating factor of e^{−ik^{3}t}:

e^{−ik^{3}t} \frac{d}{dt}û_{k} + e^{−ik^{3}t}f̂_{k} − ik^{3}e^{−ik^{3}t}û_{k} = 0

⟹ \frac{d}{dt} (e^{−ik^{3}t}û_{k}) + e^{−ik^{3}t}f̂_{k} = 0

We use the fourth order Runge-Kutta as presented above in the time step analysis to solve this problem. To analyze the stability we ignore the non-linear term uu_{x}. With this term ignored the differential equation in the Fourier domain is

û_{t} = ik^{3}û

giving us a value of λ = ik^{3} when compared to the model problem. Looking at the stability region for the forth order Runge-Kutta method we see that we require

|k^{3}Δt| ≲ 2.8 ⇒ Δt ≲ \frac{2.8}{k^{3}} ≤ \frac{2.8 · 8}{N^{3}} = \frac{22.4}{N^{3}}.

This is in good agreement with the empirical results. Figure 1(c) presents plots of the solution evolving over time for values of Δt giving stable and unstable results.

\begin{matrix}Δt = 25/N^{3} \\ Δt = 30/N^{3} \end{matrix}

Figure 1(c): u_{t} + uu_{x} + u_{xxx} = 0 with N = 256 (d) Schrödinger’s equation, Gross-Pitaevskii equation models the behavior of Bose-Einstein

con densates (BEC)

iϵu_{t} + ϵ^{2}u_{xx} + (−|u|^{2} − x^{2})u = 0 \{ \begin{matrix}ϵ = 0.3 \\ u(x, 0) = exp(−2x^{2}) \end{matrix}

We assume periodic boundary conditions for x ∈ [−π, π] and use the same spatial grid as in 1(a). Applying the Fourier transform to the discretized PDE, we get the following system of ODE’s:

\frac{d}{dt}û_{k} + iϵk^{2}û_{k} = \frac{i}{ϵ}f̂_{k} ,

where f(x, t) = (−|u(x, t)|^{2} − x^{2})u(x, t) and the partial derivatives are computed in the Fourier domain, using (1). Since this is a stiff ODE, we use an integrating factor of e^{iϵk^{2}t}:

e^{iϵk^{2}t} \frac{d}{dt}û_{k} + iϵe^{iϵk^{2}t}k^{2}û_{k} = \frac{i}{ϵ}e^{iϵk^{2}t}f̂_{k}

⟹ \frac{d}{dt} (e^{iϵk^{2}t}û_{k}) = \frac{i}{ϵ}e^{iϵk^{2}t}f̂_{k}

We use the fourth order Runge-Kutta as presented above in the time step analysis to solve this problem. For the stability analysis we again neglect the non-linear terms in order to relate the given pde to the model equation. In the Fourier domain we therefore consider

iϵû_{t} + ϵ^{2}(ik)^{2}û_{xx} = 0 ⇔ û_{t} = −ϵik^{2}û.

which gives λ = −ϵik^{2} when compared to the model problem. Therefore for stability we require

|λΔt| ≲ 2.8 ⇒ Δt ≲ \frac{2.8}{ϵk^{2}} ≤ \frac{2.8 · 4}{0.3N^{2}} ≈ \frac{38}{N^{2}} .

This seems to be in reasonable agreement with the results obtained empirically for N = 128. The following problems are all solved in a similar manner, differing only in the specific functions f(x, t).

(e) Nonlinear Schrödinger’s equation for super-fluids

iϵu_{t} + ϵ^{2}u_{xx} + (|u|^{2} − 1)u = 0 \{ \begin{matrix}ϵ = 0.3 \\ u(x, 0) = exp(−2x^{2}) \end{matrix}

As we have the same equation as before bar the non-linear term we again require

Δt ≲ \frac{38}{N^{2}} .

In practice we are able to achieve stability with a time step about 4 times (for N = 128) larger than that predicted.

(f) Nonlinear Schrödinger’s equation for nonlinear optics in saturable media

iϵu_{t} + ϵ^{2}u_{xx} + (1 − |u|^{2})|u|^{2}u = 0 \{ \begin{matrix}ϵ = 0.3 \\ u(x, 0) = exp(−2x^{2}) \end{matrix}

Again we require

Δt ≲ \frac{38}{N^{2}}

for stability of the pde with the non-linear term omitted. This is in very good agreement with the empirical results for N = 128.

(g) Landau-Lifshitz equation for ferro-magnetization

iϵu_{t} + ϵ^{2}u_{xx} = \frac{2u_{x}^{2}ū}{1 + |u|^{2}} \{ \begin{matrix}ϵ = 2 \\ u(x, 0) = exp(−2x^{2}) \end{matrix}

As this is the same as the previous equations we again have for the stability analysis

Δt ≲ \frac{38}{N^{2}}

which again agrees well with the empirical findings.

Figures 1(d), 1(e), 1(f), 1(g) present plots of the solution for these PDE’s evolving over time for values of Δt giving stable and unstable results. The plots here differ from the previous plots, since the solution are complex. Thus, the plots represent the magnitude of the solution and the colors represent the angle of the solution.

The experimental values for Δt differ from the predicted values, but the predicted values assume that all of the non-derivative terms behave linearly. This is not the case with these Schrödinger equations, and stability as a function of the time step depend on the nonlinear terms of each equation.

\begin{matrix}Δt = .5/N \\ Δt = 1/N \end{matrix}

Figure 1(d): iϵu_{t} + ϵ^{2}u_{xx}(−|u|^{2} − x^{2})u = 0 with N = 128

\begin{matrix}Δt = 12/N \\ Δt = 13/N \end{matrix}

Figure 1(e): iϵu_{t} + ϵ^{2}u_{xx} + (|u|^{2} − 1)u = 0 with N = 128

\begin{matrix}Δt = 2/N \\ Δt = 3/N \end{matrix}

Figure 1(f): iϵu_{t} + ϵ^{2}u_{xx} + (1 − |u|^{2})|u|^{2}u = 0 with N = 128

\begin{matrix}Δt = 1/N \\ Δt = 2/N \end{matrix}

Figure 1(g): iϵu_{t} + ϵ^{2}u_{xx} = \frac{2u_{x}^{2}ū}{1+|u|^{2}} with N = 128 (h) Klein-Gordon equation:

ϕ_{tt} − ϕ_{xx} + V′(ϕ) = 0  \{ \begin{matrix}V (ϕ) = \frac{1}{2}ϕ^{2} + σϕ^{4} \\ σ > 0 is a constant \\ ϕ(x, 0) = sech ((x − 2)/\sqrt{1 − c^{2}}) ,  c = .98 \\ ϕ_{t}(x, 0) = 0 \end{matrix}

We assume periodic boundary conditions for x ∈ [−π, π] and use the same spatial grid as in 1(a). Applying the Fourier transform to the discretized PDE, we get the following system of ODE’s:

\frac{d^{2}}{dt^{2}} ϕ̂_{k} + k^{2}ϕ̂_{k} + f̂_{k} = 0 ,

where f(x, t) = V′(ϕ(x, t)) and the partial derivatives are computed in the Fourier domain, using (1).

Using the second order leap frog method, we have the following iteration

\frac{ϕ̂_{k}^{n+1} − 2ϕ̂_{k}^{n} + ϕ̂_{k}^{n−1}}{(Δt)^{2}} = −f̂_{k}^{n} − k^{2}ϕ̂_{k}^{n}

with ϕ̂_{k}^{0} = F(ϕ(x, 0)), the Fourier transform of the given initial condition. To start this implicit method, use the initial condition ϕ_{t}(x, 0) = 0 and set ϕ̂_{k}^{−1} = ϕ̂_{k}^{0}, since ϕ is initially not changing with respect to t.

If we ignore the non-linear terms this pde reduces to the wave equation, therefore as we are using the 2nd order leap frog method the stability analysis is the same as that for question (1b), i.e.

Δt < \frac{4}{N}

which can be seen to be in good agreement with the empirical observations. Figure 1(h) presents plots of the solution evolving over time for values of Δt giving stable and unstable results.

\begin{matrix}Δt = 3.84/N \\ Δt = 3.85/N \end{matrix}

Figure 1(h): ϕ_{tt} − ϕ_{xx} + V′(ϕ) = 0 with N = 128 (2) Explore spectral method and nonlinear dissipation type equations.

Presented below are spectral codes to solve the following time-dependent PDEs with periodic boundary conditions in [−π, π].

As in the previous problem, we assume periodic boundary conditions for x ∈ [−π, π] and use the same spatial grid as in 1(a). The systems of ODE’s have the following form:

\frac{d}{dt}û_{k} = f̂_{k} + c(k)û_{k}

where f(x, t) is different for each equation, c(k) is dependent on the linear partial derivative terms in each equation, and the partial derivatives are computed in the Fourier domain, using (1).

We now apply the Crank-Nicolson method to solve this system:

\frac{û_{k}^{n+1} − û_{k}^{n}}{Δt} = f̂_{k}^{n+ \frac{1}{2}} + c(k)û_{k}^{n+ \frac{1}{2}} .

By expanding f̂_{k}^{n} and f̂_{k}^{n−1} in a Taylor series, we can approximate the nonlinear terms f̂_{k}^{n+ \frac{1}{2}} as follows:

f̂_{k}^{n+ \frac{1}{2}} ≈ \frac{1}{2} (3f̂_{k}^{n} − f̂_{k}^{n−1})

We then use a simple average for û_{k}^{n+ \frac{1}{2}} , since this is a linear term. Thus, the implementation of the Crank-Nicolson method that we use throughout this problem uses the following iteration:

\frac{û_{k}^{n+1} − û_{k}^{n}}{Δt} = \frac{1}{2} (3f̂_{k}^{n} − f̂_{k}^{n−1}) + c(k)\frac{û_{k}^{n+1} + û_{k}^{n}}{2} .

In question 2 the Crank-Nicholson method is used exclusively. Applying this to the model problem gives

(1 − \frac{1}{2}Δtλ)u^{n+1} − (1 + \frac{1}{2}Δtλ)u^{n} = 0

form which we infer the stability polynomial to be

π(u) = (1 − \frac{1}{2}Δtλ)u − (1 + \frac{1}{2}Δtλ).

Therefore, for absolute stability we require

|\frac{1 + \frac{1}{2}Δtλ}{1 − \frac{1}{2}Δtλ}| < 1.

Problems (a) and (c) when non-linear terms are ignored reduce to

û_{t} = −ϵk^{2}û

giving λ = −ϵk^{2}. We can see that

|\frac{1 − \frac{1}{2}ϵk^{2}Δt}{1 + \frac{1}{2}ϵk^{2}Δt}| < 1.

will be true for all Δt > 0 and therefore the method should always be stable.

Problems (d), (f), and (g) become

û_{t} = −ϵ^{2}k^{4}û

when the non-linear terms are dropped. This gives λ = −ϵ^{2}k^{4} and we see that

|\frac{1 − \frac{1}{2}ϵ^{2}k^{4}Δt}{1 + \frac{1}{2}ϵ^{2}k^{4}Δt}| < 1.

is again stable for all Δt > 0.

In problem (b) the PDE, with the non-linear terms omitted, becomes in the Fourier domain

û_{t} = k^{2}(1 − ϵ^{2}k^{2})û

so λ = k^{2}(1 − ϵ^{2}k^{2}). If the stability condition is to be met we require λ < 0. We therefore require

1 < ϵ^{2}k^{2} < \frac{ϵ^{2}N^{2}}{4}

where we have taken k_{max} = N/2. With ϵ = 0.2 this gives

N^{2} > \frac{4}{0.2^{2}}

⇔ N > 10.

Providing this condition is met we should always get stability.

As problems (h) and (g) are entirely non-linear we can not simplify them to the model equation to get a prediction for the time step required for stability. Stability analysis for problem (i) follows the analysis of the second order leap frog method in problem 1, as this problem is solved with this method. However, due to the significant nonlinear terms in (i), the experimental results do not closely reflect the predicted constraints.

(a) Burger’s equation and shock waves

u_{t} + uu_{x} = ϵu_{xx}, ϵ = 0.15, u(x, 0) = sin^{2}(x) for x < 0, u(x, 0) = 0 for x > 0 .

Δt = 1/N

Figure 2(a): u_{t} + uu_{x} = ϵu_{xx} with N = 128 (b) Kuramoto-Sivashinsky equation

u_{t} + uu_{x} = −u_{xx} − ϵ^{2}u_{xxxx}, ϵ = 0.2u(x, 0) = exp(−10x^{2}) .

Δt = 1/N

Figure 2(b): u_{t} + uu_{x} = −u_{xx} − ϵ^{2}u_{xxxx} with N = 128 (c) Allen-Cahn equation

u_{t} = u(1 − u^{2}) + ϵ^{2}u_{xx}, ϵ = 0.2, u(x, 0) = 2 sin(4x) exp(sin(x)) .

\begin{matrix}u(x, t),  Δt = 4/N \\ u_{t}(x, t),  Δt = 4/N \end{matrix}

Figure 2(c): u_{t} = u(1 − u^{2}) + ϵ^{2}u_{xx} with N = 128 (d) Cahn-Hilliard equation

u_{t} = −((1 − u^{2})u_{x})_{x}, −ϵ^{2}u_{xxxx} ϵ = 0.2, u(x, 0) = 2 sin(4x) exp(sin(x)) .

\begin{matrix}u(x, t),  Δt = 1/N \\ u_{t}(x, t),  Δt = 1/N \end{matrix}

Figure 2(d): u_{t} = −((1 − u^{2})u_{x})_{x}, −ϵ^{2}u_{xxxx} with N = 128 (f) Thin solid fluid film with slop selection

u_{t} = ∂_{x}((a(u_{x})^{2} − 1)u_{x}) − ϵ^{2}u_{xxxx}, , a = 2, ϵ = 0.2

with initial data

u_{0}(x) = 0.03 ∑_{j=1}^{N} sin(k_{j}x + c_{j}), k_{1} = 2, k_{2} = 5, k_{3} = 7,

Δt = 1/N

Figure 2(f): u_{t} = ∂_{x}((a(u_{x})^{2} − 1)u_{x}) − ϵ^{2}u_{xxxx} with N = 128 (g) Molecular Beam Expitaxy equation

u_{t} = −∂_{x}(\frac{u_{x}}{1 + (u_{x})^{2}} ) − ϵ^{2}u_{xxxx}, ϵ = 0.2, u(x, 0) = exp(−10x^{2})

Δt = 1/N

Figure 2(g): u_{t} = −∂_{x}(\frac{u_{x}}{1+(u_{x})^{2}} ) − ϵ^{2}u_{xxxx} with N = 128 (h) Thin solid film growing with slope selection

u_{t} = −ϵ^{2}(uu_{xxx})_{x} − (u^{3}u_{x})_{x}, ϵ = 0.2, u(x, 0) = exp(−x^{2})

Δt = 3.6/N

Figure 2(h): u_{t} = −ϵ^{2}(uu_{xxx})_{x} − (u^{3}u_{x})_{x} with N = 128 (i) Viscoelastic equation for phase transition:

u_{tt} − (W′(u_{x}))_{x} − u_{xxt} + u = 0, u(x, 0) = exp(−x^{2});

where W (.) is a double-well potential, say W (s) = \frac{1}{4}(s^{2} − 1)^{2}.

Δt = 1/N

Figure 2(i): u_{tt} − (W′(u_{x}))_{x} − u_{xxt} + u = 0 with N = 128 (3) Explore Chebyshev spectral method and boundary value problems.

The Chebyshev spectral method for boundary value problems differs from the Fourier spectral method in the following two ways.

• Algebraic polynomials are used instead of trigonometric polynomials.

• The grid points, rather than being equidistant, satisfy

x_{j} = cos(jπ/N), j = 0, 1, …, N.

These are referred to as Chebyshev points.

With these points in mind the Chebyshev spectral derivative can be obtained using the following procedure.

• Let p be the unique polynomial of degree N or less that satisfies p(x_{j}) = u_{j}, 0 ≤ j ≤ N, where u_{j} is an approximant to the value of the unknown function at the grid points.

• Let w_{j} = p′(x_{j}).

This last statement can be written as the linear system

w = D_{N}u

where D_{N} is the Chebyshev differentiation matrix used to approximate u′ on the discrete grid. The boundary conditions associated with the partial differential equation we wish to solve can be imposed by suitable modification of D_{N}. We also note that higher spectral derivatives can be obtained by repeatedly applying D_{N}, e.g. the second spectral derivative is given by D_{N}^{2}.

(a) We are asked to solve

u_{xx} = e^{u}, u(±1) = 0  (2)

in the interval (−1, 1). Discretizing and using the second spectral derivative gives

D̃_{N}^{2}u = exp(u)

where D̃_{N}^{2} is obtained from D_{N}^{2} by removing the the first and last rows and columns to impose the boundary conditions. Therefore we are left with a system of equations for u_{1}, …, u_{N−1} with u_{0} = u_{N} = 0. Because of the non-linear term on the right this system can not be solved by calculating the inverse of D_{N}^{2} or by simply applying the backslash operator, therefore we implement the iterative scheme

D̃_{N}^{2}u^{(n+1)} = exp(u^{(n)})

with initial guess u = 0. This method is implemented in p3a.m which corresponds to program 14 in Trefethen’s book. The output is shown in Figure 3(a).

Figure 3(a): Spectral solution of (2).

(b) In this question we are asked to compute the eigenvalues of

u_{xx} = λu, u(±1) = 0  (3)

in (−1, 1). The eigenvalues and eigenvectors are known to be

λ = −\frac{π^{2}n^{2}}{4} , n = 1, 2, …

and

u_{n} = sin (\frac{nπ(x + 1)}{2}), n = 1, 2, …

respectively. Discretizing and using the spectral derivative we get

D̃_{N}^{2}u = λu

where D̃_{N}^{2} is the same as before. We can hence use the matlab function eig to calculate the eigenvalues and eigenvectors of D̃_{N}^{2} to obtain approximations to the eigenvalues and eigenvectors of (3). Program p3b.m implements this method (corresponding to program 15 in Trefethen’s book). The results are shown in Figure 3(b). The results are accurate for λ of small amplitude but become increasingly inaccurate as λ gets larger. The reason for this is that the larger the value of the eigenvalue the more oscillatory is its corresponding eigenvector and the more oscillatory eigenvectors cannot be resolved on the discritized grid.

(c) Here we are asked to solve

u_{xx} + u_{yy} = 10 sin(8x(y − 1)), u = 0 on ∂Ω  (4)

with Ω = (0, 1) × (0, 1). We discretize by using Chebyshev points in both the x and the y directions and impose boundary conditions just as before by deleting the last rows and columns of the Chebyshev differentiation matrices. We note that second Chebyshev differentiation

Figure 3(b): Spectral solution of (3).

matrices in 2-d can be obtained using Kronecker products so the discrete system becomes

(I ⊗ D̃_{N}^{2} + D̃_{N}^{2} ⊗ I)u = f

where u is a vector containing the unknowns obtained by labeling the grid points lexigraphically and f is the discretized right hand side evaluated at the discrete grid points. This method is implemented in p3c.m corresponding to program 16 in Trefethen’s book. The result is shown in Figure 3(c).

Figure 3(c): Spectral solution to (4).

A Code Listings

% p1a.m
% 1(a) variable coefficient wave equation
% PDE: u_t + c(x)*u_x = 0
%      c(x) = .2+sin(x+2).^2, u(x,0) = exp(-100*(x+2).^2)
%      x in [-pi, pi]
% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2);               % spatial grid
t = 0;                            % initial time value
dt = 1.82/N;                      % stable time step
dt = 1.83/N;                      % unstable time step
c = .2 + sin(x+2).^2;             % c(x)
v = exp(-100*(x+2).^2);           % ~ u(x,0)
vold = exp(-100*(x-0*dt+2).^2);   % ~ u(x,-dt)
% Time-stepping by leap frog formula:
tmax = 8;                     % end time value
tplot = .15;                  % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [v; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots
    t = t+dt;                 % increment time
    v_hat = fft(v);           % differntiation via FFT
    w_hat = 1i*[0:N/2-1 0 -N/2+1:-1] .* v_hat; % discrete diff. of v
    w = real(ifft(w_hat));    % ~ u_x
    vnew = vold - 2*dt*c.*w;  % leap frog approximation
    vold = v;                 % update v_{n-1}
    v = vnew;                 % update v_n
end
  data(i+1,:) = v;            % include solution for plotting
  tdata = [tdata; t];         % time of current solution
  if max(abs(v)>2.5)          % check for instability
      data = data(1:i+1,:);   % stop when instability dominates too much
      break
  end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -3 3]), xlabel x, ylabel t, zlabel u, grid off
close(h),figure(H);

% p1b.m
% 1(b) 2nd order wave
% PDE: u_tt = u_xx
%      u(x,0) = exp(-10x^2), u_t(x,0) = 0
%      x in [-pi, pi]
% Grid, variable coefficient, and initial data:
N = 80;                          % number of grid points
h = 2*pi/N;                      % spatial step
x = h*(-N/2+1:N/2);              % spatial grid
t = 0;                           % initial time value
dt = 3.9/N;                      % stable time step
dt = 4/N;                        % unstable time step
v = exp(-10*(x.^2));             % ~ u(x,0)
vold = v;                        % ~ u(x,-dt) since u_t(x,0)=0
% Time-stepping by Leap Frog formula:
tmax = 8;                        % end time value
tplot = .075;                    % time step for plots
clf, drawnow, set(gcf,'renderer','zbuffer') % graphics rendering info
plotgap = round(tplot/dt);       % time gap between plots
dt = tplot/plotgap;              % adjust dt to line up with time plots
nplots = round(tmax/tplot);      % number of time values to plot
data = [v; zeros(nplots,N)];     % approximate solutions for plotting
tdata = t;                       % keep track of time values used
H = gcf;                         % get handle to current figure
clf, drawnow;                    % clear figure
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots                 % loop through all time steps to be plotted
    waitbar(i/nplots)            % update wait bar
    for n = 1:plotgap            % loop through all time steps between plots
        t = t+dt;                % increment time
        v_hat = fft(v);          % differntiation via FFT
        w_hat = (1i*[0:N/2-1 N/2 -N/2+1:-1]).^2 .* v_hat; % scale v_hat
        w = real(ifft(w_hat));   % ~ u_xx
        vnew = 2*v-vold+dt^2*w;  % leap frog approximation
        vold = v;                % update v_{n-1}
        v = vnew;                % update v_n
    end
    data(i+1,:) = v;             % spatial data
    tdata = [tdata; t];          % time data
    if max(abs(v)>2.5)           % check for instability
        data = data(1:i+1,:);    % save data before instability
        break
    end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -2 2]), xlabel x, ylabel t, zlabel u, grid off
close(h);figure(H);

% p1c.m
% p1c.m - Solve KdV eq. u_t + uu_x + u_xxx = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(-ik^3t)*u-hat.
%         u(x,0) = 3A^2sech^2(A(x+2)/2)+3B^2sech^2(B(x+1)/2),A=25,B=16
% Set up grid and two-soliton initial data:
N = 256;                           % number of spatial grid points
h = 2*pi/N;                        % spatial grid size
x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
dt = 25/N^3;                       % stable time step
dt = 30/N^3;                       % unstable time step
A = 25; B = 16;                    % parameters for initial value
u = 3*A^2*sech(.5*(A*(x+2))).^2 + ...
    3*B^2*sech(.5*(B*(x+1))).^2;   % initial value of u
v = fft(u);
k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
ikv = (1i*k).^3;                   % scaling factor for 3rd partial deriv.
% Solve PDE and plot results:
tmax = 0.006;                      % end time
nplt = floor((tmax/25)/dt);        % number of solutions to plot
nmax = round(tmax/dt);             % number of solutions to compute
udata = u;                         % keep track of solutions
tdata = 0;                         % keep track of time steps
clf, drawnow;                      % clear plot
set(gcf,'renderer','zbuffer');     % set rendering more for fast plotting
H = gcf;                           % get handle of current plot
h = waitbar(0,'please wait...');   % create waitbar
t = 0;
E = exp(t*ikv);                    % integrating factor at t
Edt2 = exp(-dt*ikv/2);             % integrating factor at dt/2
Edt = exp(-dt*ikv);                % integrating factor at dt
for n = 1:nmax
   t = n*dt;                       % set time for solution
   g = -1;                         % scaling factor of RHS
   % 4th-order Runge-Kutta method
   F1 = dt*g.*fft(f1c(v,k));
   F2 = dt*g.*Edt2.*fft(f1c(v+F1/2,k));
   F3 = dt*g.*Edt2.*fft(f1c(v+F2/2,k));
   F4 = dt*g.*Edt.*fft(f1c(v+F3,k));
   v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;
   if mod(n,nplt) == 0             % plot the solution
    u = real(ifft(v));             % compute solution in real domain
    waitbar(n/nmax);               % update waitbar
    udata = [udata u];             % keep track of plotted data
    tdata = [tdata t];             % keep track of plotted time
  end
end
% Trefethen's code for plotting
waterfall(x,tdata,udata'), colormap(1e-6*[1 1 1]); view(-20,25)
xlabel x, ylabel t, axis([-pi pi 0 tmax 0 2000]), grid off
set(gca,'ztick',[0 2000]), close(h), pbaspect([1 1 .13]),figure(H);

% f1c.m
function f = f1c(v,k)
v_0 = real(ifft(v));
v_x = real(ifft(i*k.*v));
f =  v_0.*v_x;

% p1d.m
% p1d.m - Gross-Pitaevskii equation for Bose-Einstein condensates
%         ieu_t + e^2u_xx + (|u|^2-x^2)u = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(iek^2t)*u_hat.
%         u(x,0) = exp(-2x^2), e = 0.3
% Set up grid and two-soliton initial data:
N = 128;                           % number of spatial grid points
h = 2*pi/N;                        % spatial grid size
x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
dt = .5/N;                         % stable time step
dt = 1/N;                          % unstable time step
u = exp(-2*(x.^2));                % initial condition
e = 0.3;                           % parameter epsilon
t = 0;                             % initial time
v = fft(u);                        % Fourier transform of u (u_hat)
k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
ikv = e*1i*k.^2;                   % scaling factor for 2nd partial deriv.
% Solve PDE and plot results:
tmax = 8;                          % end time
tplot = .1;                        % time step for plots
clf, drawnow;                      % clear figure
set(gcf,'renderer','zbuffer')      % for fast plotting
plotgap = round(tplot/dt);         % time gap between plots
dt = tplot/plotgap;                % adjust dt to line up with time plots
nplots = round(tmax/tplot);        % number of time values to plot
udata = [u zeros(N,nplots)];       % approximate solutions for plotting
tdata = t;                         % keep track of time values used
H = gcf;                           % get handle to current figure
clf, drawnow;                      % clear figure
h = waitbar(0,'please wait...');   % create waitbar
E = exp(t*ikv);                    % integrating factor at t
Edt2 = exp(dt*ikv/2);              % integrating factor at dt/2
Edt = exp(dt*ikv);                 % integrating factor at dt
for i = 1:nplots                   % loop through all time steps to be plotted
  waitbar(i/nplots)                % update wait bar
  for n = 1:plotgap                % loop through all time steps between plots
    t = t+dt;                      % set time for solution
    g = 1i/e;                     % scaling factor of RHS
    % 4th-order Runge-Kutta method
    F1 = dt*g.*fft(f1d(v,k,x));
    F2 = dt*g.*Edt2.*fft(f1d(v+F1/2,k,x));
    F3 = dt*g.*Edt2.*fft(f1d(v+F2/2,k,x));
    F4 = dt*g.*Edt.*fft(f1d(v+F3,k,x));
    v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;
  end

  u = ifft(v);                    % compute solution for plotting
  udata(:,i+1) = u;               % keep track of plotted data
  tdata = [tdata t];              % keep track of plotted time
  if max(abs(u)>5)                % check for instability
     udata = udata(:,1:i+1);      % save data before instability
     break
 end
end
% Trefethen's code for plotting
waterfall(x,tdata,abs(udata'),angle(udata')); % for complex values
view(-20,25), xlabel x, ylabel t, zlabel u,
axis([-pi pi 0 tmax min(min(udata)) 1]), grid off
set(gca,'ztick',[-1 1]), close(h), pbaspect([1 1 .5]),figure(H);
% f1d.m
function f = f1d(v,k,x)
v_0 = ifft(v);
f =  (-abs(v_0).^2-x.^2).*v_0;

% p1e.m
% p1e.m - Super fluids
%         ieu_t + e^2u_xx + (|u|^2-1)u = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(iek^2t)*u_hat.
%         u(x,0) = exp(-2x^2), e = 0.3
% Set up grid and two-soliton initial data:
N = 128;                           % number of spatial grid points
h = 2*pi/N;                        % spatial grid size
x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
dt = 12/N;                         % stable time step
dt = 13/N;                         % unstable time step
u = exp(-2*(x.^2));                % initial condition
e = 0.3;                           % parameter epsilon
t = 0;                             % initial time
v = fft(u);                        % Fourier transform of u (u_hat)
k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
ikv = e*1i*k.^2;                   % scaling factor for 2nd partial deriv.
% Solve PDE and plot results:
tmax = 8;                         % end time
tplot = .15;                        % time step for plots
clf, drawnow;                      % clear figure
set(gcf,'renderer','zbuffer')      % for fast plotting
plotgap = round(tplot/dt);         % time gap between plots
dt = tplot/plotgap;                % adjust dt to line up with time plots
nplots = round(tmax/tplot);        % number of time values to plot
udata = [u zeros(N,nplots)];       % approximate solutions for plotting
tdata = t;                         % keep track of time values used
H = gcf;                           % get handle to current figure
clf, drawnow;                      % clear figure
h = waitbar(0,'please wait...');   % create waitbar
E = exp(t*ikv);                    % integrating factor at t

Edt2 = exp(dt*ikv/2);              % integrating factor at dt/2
Edt = exp(dt*ikv);                 % integrating factor at dt
for i = 1:nplots                   % loop through all time steps to be plotted
  waitbar(i/nplots)                % update wait bar
  for n = 1:plotgap                % loop through all time steps between plots
    t = t+dt;                      % set time for solution
    g = 1i/e;                     % scaling factor of RHS
    % 4th-order Runge-Kutta method
    F1 = dt*g.*fft(f1e(v,k));
    F2 = dt*g.*Edt2.*fft(f1e(v+F1/2,k));
    F3 = dt*g.*Edt2.*fft(f1e(v+F2/2,k));
    F4 = dt*g.*Edt.*fft(f1e(v+F3,k));
    v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;
  end
  u = ifft(v);                    % compute solution for plotting
  udata(:,i+1) = u;               % keep track of plotted data
  tdata = [tdata t];              % keep track of plotted time
  if max(abs(u)>2.5)              % check for instability
     udata = udata(:,1:i+1);      % save data before instability
     break
  end
end
% Trefethen's code for plotting
waterfall(x,tdata,abs(udata'),angle(udata')); % for complex values
view(-20,25), xlabel x, ylabel t, zlabel u, axis([-pi pi 0 tmax 0 1]);
grid off, set(gca,'ztick',[-1 1]), close(h), pbaspect([1 1 .5]),figure(H);
% f1e.m
function f = f1e(v,k)
v = ifft(v);
f =  (abs(v).^2-1).*v;

% p1f.m
% p1f.m - Nonlinear optics
%         ieu_t + e^2u_xx + (1-|u|^2)|u|^2u = 0 on [-pi,pi] by
%         FFT with integrating factor v = exp(iek^2t)*u_hat.
%         u(x,0) = exp(-2x^2), e = 0.3
% Set up grid and two-soliton initial data:
N = 128;                           % number of spatial grid points
h = 2*pi/N;                        % spatial grid size
x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
dt = 2/N;                          % stable time step
dt = 3/N;                          % unstable time step
u = exp(-2*(x.^2));                % initial condition
e = 0.3;                           % parameter epsilon
t = 0;                             % initial time
v = fft(u);                        % Fourier transform of u (u_hat)
k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
ikv = e*1i*k.^2;                   % scaling factor for 2nd partial deriv.

% Solve PDE and plot results:
tmax = 8;                          % end time
tplot = .15;                        % time step for plots
clf, drawnow;                      % clear figure
set(gcf,'renderer','zbuffer')      % for fast plotting
plotgap = round(tplot/dt);         % time gap between plots
dt = tplot/plotgap;                % adjust dt to line up with time plots
nplots = round(tmax/tplot);        % number of time values to plot
udata = [u zeros(N,nplots)];       % approximate solutions for plotting
tdata = t;                         % keep track of time values used
H = gcf;                           % get handle to current figure
clf, drawnow;                      % clear figure
h = waitbar(0,'please wait...');   % create waitbar
E = exp(t*ikv);                    % integrating factor at t
Edt2 = exp(dt*ikv/2);              % integrating factor at dt/2
Edt = exp(dt*ikv);                 % integrating factor at dt
for i = 1:nplots                   % loop through all time steps to be plotted
  waitbar(i/nplots)                % update wait bar
  for n = 1:plotgap                % loop through all time steps between plots
    t = t+dt;                      % set time for solution
    g = 1i/e;                     % scaling factor of RHS
    % 4th-order Runge-Kutta method
    F1 = dt*g.*fft(f1f(v,k));
    F2 = dt*g.*Edt2.*fft(f1f(v+F1/2,k));
    F3 = dt*g.*Edt2.*fft(f1f(v+F2/2,k));
    F4 = dt*g.*Edt.*fft(f1f(v+F3,k));
    v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;
  end
  u = ifft(v);                    % compute solution for plotting
  udata(:,i+1) = u;               % keep track of plotted data
  tdata = [tdata t];              % keep track of plotted time
  if max(abs(u)>2.5)              % check for instability
     udata = udata(:,1:i+1);      % save data before instability
     break
  end
end
% Trefethen's code for plotting
waterfall(x,tdata,abs(udata'),angle(udata'));  %for complex values
view(-20,25), xlabel x, ylabel t, zlabel u, axis([-pi pi 0 tmax 0 1])
grid off, set(gca,'ztick',[-1 1]), close(h), pbaspect([1 1 .5]),figure(H);
% f1f.m
function f = f1f(v,k)
v = ifft(v);
f =  (1-abs(v).^2).*(abs(v).*2).*v;

% p1g.m
% p1g.m - Landau - Lifshitz equation for ferro-magnetization
%         ieu_t + e^2u_xx = 2*(u_x)^2*uconj/(1+|u|^2) on [-pi,pi] by
%         FFT with integrating factor v = exp(iek^2t)*u_hat.
%         u(x,0) = exp(-2x^2), e = 2

% Set up grid and two-soliton initial data:
N = 128;                           % number of spatial grid points
h = 2*pi/N;                        % spatial grid size
x = (2*pi/N)*(-N/2+1:N/2)';        % spatial grid
dt = 1/N;                          % stable time step
dt = 2/N;                          % unstable time step
u = exp(-2*(x.^2));                % initial condition
e = 2;                             % parameter epsilon
t = 0;                             % initial time
v = fft(u);                        % Fourier transform of u (u_hat)
k = [0:N/2-1 0 -N/2+1:-1]';        % permuted k to match freq. of fft
ikv = e*1i*k.^2;                   % scaling factor for 3rd partial deriv.
% Solve PDE and plot results:
tmax = 3;                         % end time
tplot = .05;                        % time step for plots
clf, drawnow;                      % clear figure
set(gcf,'renderer','zbuffer')      % for fast plotting
plotgap = round(tplot/dt);         % time gap between plots
dt = tplot/plotgap;                % adjust dt to line up with time plots
nplots = round(tmax/tplot);        % number of time values to plot
udata = [u zeros(N,nplots)];       % approximate solutions for plotting
tdata = t;                         % keep track of time values used
H = gcf;                           % get handle to current figure
clf, drawnow;                      % clear figure
h = waitbar(0,'please wait...');   % create waitbar
E = exp(t*ikv);                    % integrating factor at t
Edt2 = exp(dt*ikv/2);              % integrating factor at dt/2
Edt = exp(dt*ikv);                 % integrating factor at dt
for i = 1:nplots                   % loop through all time steps to be plotted
  waitbar(i/nplots)                % update wait bar
  for n = 1:plotgap                % loop through all time steps between plots
    t = t+dt;                      % set time for solution
    g = -1i/e;                     % scaling factor of RHS
    % 4th-order Runge-Kutta method
    F1 = dt*g.*fft(f1g(v,k));
    F2 = dt*g.*Edt2.*fft(f1g(v+F1/2,k));
    F3 = dt*g.*Edt2.*fft(f1g(v+F2/2,k));
    F4 = dt*g.*Edt.*fft(f1g(v+F3,k));
    v = Edt.*v + (Edt.*F1 + 2*Edt2.*(F2+F3) + F4)/6;
  end
  u = ifft(v);              % compute solution for plotting
  udata(:,i+1) = u;               % keep track of plotted data
  tdata = [tdata t];              % keep track of plotted time
  if max(abs(u)>5)              % check for instability
     udata = udata(:,1:i+1);      % save data before instability
     break
  end
end
% Trefethen's code for plotting
waterfall(x,tdata,abs(udata'),angle(udata')); view(-20,25);
xlabel x, ylabel t, zlabel u, axis([-pi pi 0 tmax min(min(udata)) 2]);
grid off, set(gca,'ztick',[-1 1]), close(h), pbaspect([1 1 .5]),figure(H);

% f1g.m
function f = f1g(v,k)
v_0 = ifft(v);
v_x = ifft(i*k.*v);
v_c = ifft(conj(fliplr(v)));
f =  2*(v_x.^2).*v_c./(1+abs(v_0).^2);

% p1h.m
% 1(h) Klein-Gordon
% PDE: phi_tt - phi_xx + V'(phi) , V(phi) = .5*phi^2+s*phi^4
%      phi(x,0) = sech((x-2)/sqrt(1-c^2)), phi_t(x,0) = 0
%      x in [-pi, pi]
% Grid, variable coefficient, and initial data:
N = 128;                         % number of grid points
h = 2*pi/N;                      % spatial step
x = h*(-N/2+1:N/2);              % spatial grid
t = 0;                           % initial time value
dt = 3.84/N;                     % stable time step
dt = 3.85/N;                     % unstable time step
s = 10;                          % sigma for V'(phi)
c = .98;                         % parameter of initial condition
v = sech((x-2)/sqrt(1-c^2));     % initial condition
vold = v;                        % ~ u(x,-dt) since u_t(x,0)=0
% Time-stepping by Leap Frog formula:
tmax = 20;                       % end time value
tplot = .075;                    % time step for plots
clf, drawnow, set(gcf,'renderer','zbuffer') % graphics rendering info
plotgap = round(tplot/dt);       % time gap between plots
dt = tplot/plotgap;              % adjust dt to line up with time plots
nplots = round(tmax/tplot);      % number of time values to plot
data = [v; zeros(nplots,N)];     % approximate solutions for plotting
tdata = t;                       % keep track of time values used
H = gcf;                         % get handle to current figure
clf, drawnow;                    % clear figure
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots                 % loop through all time steps to be plotted
    waitbar(i/nplots)            % update wait bar
    for n = 1:plotgap            % loop through all time steps between plots
        t = t+dt;                % increment time
        v_hat = fft(v);          % differntiation via FFT
        w_hat = (1i*[0:N/2-1 N/2 -N/2+1:-1]).^2 .* v_hat; % scale v_hat
        w = real(ifft(w_hat));   % ~ u_xx
        f = v + 4*s*v.^3;
        vnew = 2*v-vold+dt^2*(w-f); % leap frog approximation
        vold = v;                % update v_{n-1}
        v = vnew;                % update v_n
    end
    data(i+1,:) = v;             % spatial data
    tdata = [tdata; t];          % time data
    if max(abs(v)>2.5)           % check for instability
        data = data(1:i+1,:);    % save data before instability

        break
    end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(-20,25),colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax min(min(data)) 1]), pbaspect([1 1 .5]),
xlabel x, ylabel t, zlabel u, grid off, close(h);figure(H);

% p2a.m
% 2(a) Burger's equation and shock waves
% PDE: u_t + u*u_x = -e*u_xx
%  IC: u(x,0) = /sin(x)^2 for x < 0
%               \   0     for x > 0
%      x in [-pi, pi]
% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';              % spatial grid
t = 0;                            % initial time value
dt = 1/N;                         % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = .15;
% initial condition: define u = u(x,0) and uold = u(x,-dt)
for i = 1:N
    if x(i) < 0
      u(i) = (sin(x(i)))^2;
      uold(i) =(sin(x(i)-.2*dt))^2;
    else
      u(i) = 0;
      uold(i) = 0;
    end
end
u_hat = fft(uold);            % step back in time
u_x =real(ifft(1i*k.*u_hat));
f_k = fft(uold.*u_x);
% Time-stepping by Crank-Nicolson
tmax = 8;                     % end time value
tplot = .1;                   % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots

t = t+dt;                 % increment time
% start creating f_k^n which I call f_k and f_k^(n-1), I call fold_k
% where f = u*u_x
uold_hat = u_hat;         % u_hat from last iteration
u_hat = fft(u);           % new u_hat
uold_x = u_x;             % u_x from last iteration
u_x =real(ifft(1i*k.*u_hat)); % new u_x
fold_k = f_k;             % f_k from last iteration
f_k = fft(u.*u_x);        % new f_k
% 2nd-order approximation of f(n+.5)
fhalf = (3*f_k - fold_k)/2;
% Now I apply Crank-Nicolson
l1 = dt^(-1)-e*k.^2/2;
l2 = dt^(-1)+e*k.^2/2;
unew_hat = (l1.*u_hat-fhalf)./l2;
unew = real(ifft(unew_hat));
uold = u;                 % update u_(n)
u = unew;                 % update u_(n+1)
end
data(i+1,:) = u;            % include solution for plotting
tdata = [tdata; t];         % time of current solution
if max(abs(u)>30)           % check for instability
  data = data(1:i+1,:);   % stop when instability dominates too much
  break
end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax 0 1]), xlabel x, ylabel t, zlabel u, grid off
close(h),figure(H);

% p2b.m
% 2(b) Kuramoto-Sivhsinsky
% PDE: u_t + u*u_x = -u_xx-e^2*u_xxxx
%  IC: u(x,0) = exp(-10*x^2)
%      x in [-pi, pi] and e = 0.2
% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';              % spatial grid
t = 0;                            % initial time value
dt = 1/N;                         % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = .2;
% initial condition: define u = u(x,0) and uold = u(x,-dt)
u = (exp(-10*x.^2))';
uold = (exp(-10*(x-.2*dt).^2))';
% Time-stepping by Crank-Nicolson

tmax = 5;                     % end time value
tplot = .05;                  % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
waitbar(i/nplots)           % update wait bar
for n = 1:plotgap           % loop through all time steps between plots
t = t+dt;                 % increment time
% start creating f_k^n which I call f_k and f_k^(n-1), I call fold_k
% where f = u*u_x
uold_hat = u_hat;         % u_hat from last iteration
u_hat = fft(u);           % new u_hat
uold_x = u_x;             % u_x from last iteration
u_x =real(ifft(1i*k.*u_hat)); % new u_x
fold_k = f_k;             % f_k from last iteration
f_k = fft(u.*u_x);        % new f_k
% 2nd-order approximation of f(n+.5)
fhalf = (3*f_k - fold_k)/2;
% Now I apply Crank-Nicolson
l1 = dt^(-1)+(k.^2-e^2*k.^4)/2;
l2 = dt^(-1)-(k.^2-e^2*k.^4)/2;
unew_hat = (l1.*u_hat-fhalf)./l2;
unew = real(ifft(unew_hat));
uold = u;                 % update v_{n-1}
u = unew;                 % update v_n
end
data(i+1,:) = u;            % include solution for plotting
tdata = [tdata; t];         % time of current solution
if max(abs(u)>30)           % check for instability
  data = data(1:i+1,:);   % stop when instability dominates too much
   break
end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax min(min(data)) 30]), xlabel x, ylabel t, zlabel u;
grid off,close(h),figure(H);

% p2c.m
% 2(c) Allen - Cahn
% PDE: u_t = u(1-u^2)+e^2*u_xx
%  IC: u(x,0) = 2*sin(4x)*exp(sin(x))

%      x in [-pi, pi] and e = 0.2
% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';              % spatial grid
t = 0;                            % initial time value
dt = 4/N;                         % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = .2;
% initial condition: define u = u(x,0) and uold = u(x,-dt)
u = (2*sin(4*x).*exp(sin(x)))';
uold = (2*sin(4*(x-.2*dt)).*exp(sin(x-.2*dt)))';
u_hat = fft(uold);
f_k = fft(uold.*(1-uold.^2));
% Time-stepping by Crank-Nicolson
tmax = 100;                   % end time value
tplot = 2;                    % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
data2 = [zeros(nplots-1,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots
    t = t+dt;                 % increment time
    % start creating f_k^n which I call f_k and f_k^(n-1), I call fold_k
    % where f = u*u_x
    uold_hat = u_hat;         % u_hat from last iteration
    u_hat = fft(u);           % new u_hat
    fold_k = f_k;             % f_k from last iteration
    f_k = fft(u.*(1-u.^2));   % new f_k
    % 2nd-order approximation of f(n+.5)
    fhalf = (3*f_k - fold_k)/2;
    % Now I apply Crank-Nicolson
    l1 = (1/dt - e^2*k.^2/2);
    l2 = (1/dt + e^2*k.^2/2);
    unew_hat = (fhalf +l1.*u_hat)./l2;
    unew = real(ifft(unew_hat));
    uold = u;                 % update v_{n-1}
    u = unew;                 % update v_n
  end
  u_xx = real(ifft(-k.^2.*unew_hat));
  uprime = unew.*(1-unew.^2)+e^2*u_xx;
  data(i+1,:) = u;            % include solution for plotting
  data2(i,:) = uprime;

  tdata = [tdata; t];         % time of current solution
  if max(abs(u)>30)           % check for instability
      data = data(1:i+1,:);   % stop when instability dominates too much
      break
  end
end
% Trefethen's plotting code
figure(1)
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -5 3]), xlabel x, ylabel t, zlabel u, grid off
close(h)
figure(2);
waterfall(x,tdata(2:end,:),data2), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -.5 .5]), xlabel x, ylabel t, zlabel u, grid off

% p2d.m
% 2(d) Cahn - Hilliard equation
% PDE: u_t = -((1-u^2)*u_x)_x - e^2*u_xxxx
%  IC: u(x,0) = 2*sin(4x)*exp(sin(x))
%      x in [-pi, pi] and e = 0.2
% Grid, variable coefficient, and initial data:
N = 128;                       % number of grid points
h = 2*pi/N;                    % spatial step
x = h*(-N/2+1:N/2)';           % spatial grid
t = 0;                         % initial time value
dt = 1/N;                      % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = 5;

% initial condition: define u = u(x,0) and uold = u(x,-dt)
u = (2*sin(4*x).*exp(sin(x)))';
uold = (2*sin(4*(x-.2*dt)).*exp(sin(x-.2*dt)))';
% Time-stepping by 2nd-order Adams-Bashford formula:
tmax = 50;                    % end time value
tplot = .5;                   % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots
    t = t+dt;                 % increment time

% start creating f_k^n which I call f_k and f_k^(n-1), I call fold_k
u_hat = fft(u);
u_x = real(ifft(i*k.*u_hat));
uaux = fft((1-u.^2).*u_x);
f_k = i*k.*uaux;
uold_hat = fft(uold);
uold_x = real(ifft(i*k.*uold_hat));
uoldaux = fft((1-uold.^2).*uold_x);
fold_k = i*k.*uoldaux;
% AB - CN
fhalf = (3*f_k - fold_k)/2;
l1 = (1/dt - e^2*k.^4/2);
l2 = (1/dt + e^2*k.^4/2);
unew_hat = (-fhalf +l1.*u_hat)./l2;
unew = real(ifft(unew_hat));
uold = u;                 % update v_{n-1}
u = unew;                 % update v_n
end
uprime = real(ifft(-f2d(u,k)))-e^2.*real(ifft(k.^4.*unew_hat));
data(i+1,:) = u;            % include solution for plotting
data2(i+1,:) = uprime;
tdata = [tdata; t];         % time of current solution
if max(abs(u)>30)          % check for instability
  data = data(1:i+1,:);   % stop when instability dominates too much
  break
end
end
% Trefethen's plotting code
figure(1)
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -6 6]), xlabel x, ylabel t, zlabel u, grid off
close(h);
figure(2);
waterfall(x,tdata,data2), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -20 300]), xlabel x, ylabel t, zlabel u, grid off
% f2d.m
function f = f2d(u,k)
u_hat = fft(u);
u_x = real(ifft(i*k.*u_hat));
uaux = fft((1-u.^2).*u_x);
f = i*k.*uaux;

% p2f.m
% 2(f) Thin solid fluid film with slop selection
% PDE: u_t = ((a(u_x)^2-1)*u_x)_x - e^2*u_xxxx
%  IC: u(x,0) = 0.03 *sum(sin(k*x+c))
%      Choose c = k, k = [2 5 7]
%      x in [-pi, pi] and e = 0.2

% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';              % spatial grid
t = 0;                            % initial time value
dt = 1/N;                         % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = 8;
a = .2;
coef= [2 5 7];
coef2 = [.1 2 1.3];
u = zeros(1,128);
uold = zeros(1,128);
% initial condition: define u = u(x,0) and uold = u(x,-dt)
for i = 1:3
      u = u + .3*sin(coef(i)*x+coef2(i))';
      uold = uold + .3*sin(coef(i)*(x-dt)+coef2(i))';
end
% Time-stepping by Crank-Nicholson
tmax = 8;                     % end time value
tplot = .1;                   % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots
    t = t+dt;                 % increment time
    % start creating f_k^n which I call f_k and f_k^(n-1), I call fold_k
    u_hat = fft(u);
    u_x = real(ifft(1i*k.*u_hat));
    u_xx = real(ifft(-k.^2.*u_hat));
    f_k = fft(3*a*u_x.^2.*u_xx-u_xx);
    uold_hat = fft(uold);
    uold_x = real(ifft(1i*k.*uold_hat));
    uold_xx = real(ifft(-k.^2.*uold_hat));
    fold_k = fft(3*a*uold_x.^2.*uold_xx-uold_xx);
    % AB - CN
    fhalf = (3*f_k - fold_k)/2;
    l1 = (1/dt - e^2*k.^4/2);
    l2 = (1/dt + e^2*k.^4/2);
    unew_hat = -(fhalf +l1.*u_hat)./l2;
    unew = real(ifft(unew_hat));
    uold = u;                 % update v_{n-1}
    u = unew;                 % update v_n
  end

  data(i+1,:) = u;            % include solution for plotting
  tdata = [tdata; t];         % time of current solution
  if max(abs(u)>2.5)          % check for instability
      data = data(1:i+1,:);   % stop when instability dominates too much
      break
  end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -1 1]), xlabel x, ylabel t, zlabel u, grid off
close(h),figure(H);

% p2g.m
% 2(g) Molecular Beam Expitaxy equation
% PDE: u_t = ((a(u_x)^2-1)*u_x)_x - e^2*u_xxxx
%  IC: u(x,0) = exp(-10x^2)
%      x in [-pi, pi] and e = 0.2
% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';                   % spatial grid
t = 0;                            % initial time value
dt = 1/N;                      % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = .2;
% initial condition: define u = u(x,0) and uold = u(x,-dt)
u = exp(-10*x.^2)';
uold = exp(-10*(x-.2*dt).^2)';
% Time-stepping by Crank-Nicolson
tmax = 2;                     % end time value
tplot = .03;                  % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots
    t = t+dt;                 % increment time
    % start creating f_k^n which I call f_k and f_k^(n-1), I call fold_k
    u_hat = fft(u);
    u_x = real(ifft(1i*k.*u_hat));
    uaux = fft(u_x./(1+(u_x).^2));

    f_k = 1i*k.*uaux;
    uold_hat = fft(uold);
    uold_x = real(ifft(1i*k.*uold_hat));
    uoldaux = fft(uold_x./(1+(uold_x).^2));
    fold_k = 1i*k.*uoldaux;
    % AB - CN
    fhalf = (3*f_k - fold_k)/2;
    l1 = (1/dt - e^2*k.^4/2);
    l2 = (1/dt + e^2*k.^4/2);
    unew_hat = (fhalf +l1.*u_hat)./l2;
    unew = real(ifft(unew_hat));
    uold = u;                 % update v_{n-1}
    u = unew;                 % update v_n
  end
  data(i+1,:) = u;            % include solution for plotting
  tdata = [tdata; t];         % time of current solution
  if max(abs(u)>2.5)          % check for instability
      data = data(1:i+1,:);   % stop when instability dominates too much
      I = find(data<-3);      % find unstable negative data
      data(I) = -3;           % set cutoff for plotting data
                              % this should = zmin in axis below
      break
  end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax -.5 .5]), xlabel x, ylabel t, zlabel u, grid off
close(h),figure(H);

% p2h.m
% 2(h) Thin solid film growing with slope selection
% PDE: u_t = -e^2*(u*u_xxx)_x-(u^3*u_x)_x
%  IC: u(x,0) = exp(-10x^2)
%      x in [-pi, pi] and e = 0.2
% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';              % spatial grid
t = 0;                            % initial time value
dt = 1;                           % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = .2;
% initial condition: define u = u(x,0) and uold = u(x,-dt)
u = exp(-x.^2)';
uold = exp(-(x-dt).^2)';
% Time-stepping by 2nd-order Adams-Bashford formula:
tmax = 500;                   % end time value
tplot = 3;                    % time step for plots

clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
for i = 1:nplots              % loop through all time steps to be plotted
waitbar(i/nplots)           % update wait bar
counter = 0;
for n = 1:plotgap           % loop through all time steps between plots
t = t+dt;                 % increment time
u_k = fft(u);
uold_k =fft(uold);
u_x = real(ifft(i*k.*u_k));
uold_x = real(ifft(i*k.*uold_k));
u_xxx = real(ifft(-1i*k.^3.*u_k));
uold_xxx = real(ifft(-1i*k.^3.*uold_k));
f1 = i*k.*fft(u.*u_xxx);
f2 = i*k.*fft(u.^3.*u_x);
f_k = -e^2*f1 - f2;
fold1 = i*k.*fft(uold.*uold_xxx);
fold2 = i*k.*fft(uold.^3.*uold_x);
fold_k = -e^2*fold1 - fold2;
fhalf = (3*f_k - fold_k)/2;
unew = uold_k - 2*dt*fhalf;
unew = real(ifft(unew));
uold = u;
u = unew;
end
data(i+1,:) = u;            % include solution for plotting
tdata = [tdata; t];         % time of current solution
if max(abs(u)>2000)         % check for instability
  data = data(1:i+1,:);   % stop when instability dominates too much
  I = find(data<-3);      % find unstable negative data
  data(I) = -3;           % set cutoff for plotting data
                          % this should = zmin in axis below
  break
end
end
% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax 0 1]), xlabel x, ylabel t, zlabel u, grid off
close(h),figure(H);

% p2i.m
% 2(i) Thin solid film growing with slope selection
% PDE: u_t = -e^2*(u*u_xxx)_x-(u^3*u_x)_x
%  IC: u(x,0) = exp(-10x^2)

% Grid, variable coefficient, and initial data:
N = 128;                          % number of grid points
h = 2*pi/N;                       % spatial step
x = h*(-N/2+1:N/2)';                   % spatial grid
t = 0;                            % initial time value
dt = 1/N;                         % stable time step
k = [0:N/2-1 0 -N/2+1:-1];
e = 5;
% initial condition: define u = u(x,0) and uold = u(x,-dt)
u = exp(-x.^2)';
uold = exp(-(x-dt).^2)';

% Time-stepping by Crank-Nicolson
tmax = 5;                     % end time value
tplot = .05;                  % time step for plots
clf, drawnow;                 % clear figure
set(gcf,'renderer','zbuffer');% set rendering for fast plotting
H = gcf;                      % get handle to current figure
plotgap = round(tplot/dt);    % time gap between plots
dt = tplot/plotgap;           % adjust dt to line up with time plots
nplots = round(tmax/tplot);   % number of time values to plot
data = [u; zeros(nplots,N)];  % approximate solutions for plotting
tdata = t;                    % keep track of time values used
h = waitbar(0,'please wait...'); % create waitbar
counter = 0;
for i = 1:nplots              % loop through all time steps to be plotted
  waitbar(i/nplots)           % update wait bar
  for n = 1:plotgap           % loop through all time steps between plots
    t = t+dt;                 % increment time
    counter = counter + 1;
    u_k = fft(u);
    uold_k =fft(uold);
    u_x = real(ifft(i*k.*u_k));
    uold_x = real(ifft(i*k.*uold_k));
    f_k = i*k.*(fft(u_x.^3-u_x));
    fold_k = i*k.*(fft(uold_x.^3-uold_x));
    fhalf = (3*f_k - fold_k)/2;
    l1 = (2/dt^2+k.^2/dt-1/2);
    l2 = (1/dt^2+k.^2/dt+1/2);
    unew = (fhalf + l1.*u_k- uold_k/dt^2)./l2;
    unew = real(ifft(unew));
    uold = u;
    u = unew;
end
  data(i+1,:) = u;            % include solution for plotting
  tdata = [tdata; t];         % time of current solution
  if max(abs(u)>100)          % check for instability
      data = data(1:i+1,:);   % stop when instability dominates too much
      break
  end
end

% Trefethen's plotting code
waterfall(x,tdata,data), view(10,70), colormap(1e-6*[1 1 1]);
axis([-pi pi 0 tmax min(min(data)).5]), xlabel x, ylabel t, zlabel u;
grid off, close(h),figure(H);

% p3a.m
% 3(a) Poisson Boltzmann equation.
% PDE: u_xx = exp(u)
%      u(-1) = u(1) = 0
%      x in [-1, 1]
% Set up Chebyshev differentiatin matrix.
N = 16;           % Number of grid points
[D,x] = cheb(N);  % Chebyshev diff. matrix
D2 = D^2;         % Second order Chebyshev diff. matrix
D2 = D2(2:N,2:N); % Remove first and last col's and rows
% Carry out iterative procedure.
u = zeros(N-1,1);             % Initial guess
change = 1;                   % Used for stopping criteria
it = 0;                       % Initialize iterative counter
while change > 1e-15          % fixed-point iteration
  unew = D2\exp(u);           % Update u
  change = norm(unew-u,inf);  % Update change
  u = unew;                   % Set old u to new u
  it = it+1;                  % Increase iterative counter
end
% Impose boundary conditions and plot solution.
u = [0;u;0];                      % Add bdy values u_0 = u_N = 0
clf                               % Clear figure
subplot('position',[.1 .4 .8 .5]) % Position of graph in fig.
plot(x,u,'.','markersize',16)     % Plot solution
xx = -1:.01:1;                    % Set up finer x grid
uu = polyval(polyfit(x,u,N),xx);  % Interp. sol'n with poly of deg. N
line(xx,uu)                       % Plot poly. interpolant
grid on                           % Put grid on graph.
title(sprintf('no. steps = %d      u(0) =%18.14f',it,u(N/2+1))) % Title

% p3b.m
% 3(b) Eigenvalue problem.
% PDE: u_xx = lambda u
%      u(-1) = u(1) = 0
%      x in [-1, 1]
% Set up Chebyshev differentiation matrix.
N = 36;                   % Number of grid points
[D,x] = cheb(N);          % Chebyshev diff. matrix
D2 = D^2;                 % Second order Chebyshev diff. matrix
D2 = D2(2:N,2:N);         % Remove first and last rows and col's

% Calulate eigenvalues.
[V,Lam] = eig(D2);        % Calulate eigenvectors and eigenvalue matrix
lam = diag(Lam);          % Eigenvalues are diagonal enteries of Lam
[foo,ii] = sort(-lam);    % Sort eigenvalues from smallest to largest
lam = lam(ii);            % Index eigenvalues
V = V(:,ii);              % Index eigenvectors
% Plot 6 eigenvectors.
clf
for j = 5:5:30
  u = [0;V(:,j);0];               % Impose bdy conditions
  subplot(7,1,j/5)                % Create subplot
  plot(x,u,'.','markersize',12),  % Plot eigenvector
  grid on                         % Put grid on graph
  xx = -1:.01:1;                  % Create finer grid
  uu = polyval(polyfit(x,u,N),xx);% Interplote soln on finer grid
  line(xx,uu),                    % Plot interpolant
  axis off                        % Remove axis' and then label plot
  text(-.4,.5,sprintf('eig %d =%20.13f*4/pi^2',j,lam(j)*4/pi^2))
  text(.7,.5,sprintf('%4.1f  ppw', 4*N/(pi*j)))
end

% p3c.m
% 3(c) Poisson equation.
% PDE: u_xx + u_yy = 10sin(8x(y-1))
%      u = 0 on bdy
%      in (0,1)x(0,1)
% Set up grids and tensor product Laplacian and solve for u:
N = 24;                             % Number of grid points in 1d
[D,x] = cheb(N);                    % Cheb. diff. matrix.
y = x;                              % Set y as Cheb. points
[xx,yy] = meshgrid(x(2:N),y(2:N));  % Obtain discrete grid
xx = xx(:);                         % Stretch 2D grids to 1D vectors
yy = yy(:);                         %              "
f = 10*sin(8*xx.*(yy-1));           % Set up RHS
D2 = D^2;                           % Second order Cheb. diff. matrix
D2 = D2(2:N,2:N);                   % Remove first and last col's
I = eye(N-1);                       % Create identity matrix
L = kron(I,D2) + kron(D2,I);        % Laplacian using kronecker product
u = L\f;                            % Solve problem
% Reshape long 1D results onto 2D grid:
uu = zeros(N+1,N+1);                % Boundary values
uu(2:N,2:N) = reshape(u,N-1,N-1);   % Insert solution and reshape
[xx,yy] = meshgrid(x,y);            % Set up discrete grid
% Interpolate to finer grid and plot:
[xxx,yyy] = meshgrid(-1:.04:1,-1:.04:1);   % Set up finer grid
uuu = interp2(xx,yy,uu,xxx,yyy,'cubic');   % Interpolate solution
mesh(xxx,yyy,uuu)                          % Plot interpolant
colormap(1e-6*[1 1 1]);                    % Color mesh lines black
xlabel x, ylabel y, zlabel u               % Label graph