Duffing Equation

Clear[ω0, ω,δ, β,γ,ϕ,x,t]
DuffEq= x''[t] + δ x'[t] + ((ω0)^2) x[t] + β ((x[t])^3)== γ Cos[ω t + ϕ]
TMax=234;
{ω0, ω,δ, β,γ,ϕ}={1,1.1, 0.02,12,4,0}
xSol[t_]=x[t]/.NDSolve[{
DuffEq,
x[0]==1.2,
x'[0]==2.1
},
x[t],
{t, 0, TMax},
MaxSteps->(10^5)][[1]];
Plot[xSol[t],{t, 0, TMax}]

KdV Equation

The idea is to work out how to implement a spectral code for a specific equation.  We are going
to do the KdV equation
	(Subscript[u, t])+(Subscript[u, xxx])- 6 u (Subscript[u, x])=0
on a periodic spatial domain.  We are going to transform in "x" and consider the spectral components.

Here is the built in Mathematica Solver. So that we can see what the solution should look like.

L= 200;
f[x_]=0.2(E^(-0.1 (((x-L/2))^2)));
TMax=23;
Eq=D[u[x,t],t]- 6 u[x,t] D[u[x,t],x] + D[u[x,t],{x,3}]==0;
TraditionalForm[Eq]
uSol=u/.NDSolve[{
Eq,
u[0,t]==u[L,t],
u[x,0]==f[x]
},
u,
{x, 0, L},
{t,0,TMax}
][[1]];
Plot3D[ uSol[x,t],
{x,0,L},
{t, 0, TMax},
PlotRange->All,
PlotPoints->100]
Manipulate[
Plot[ uSol[x,t],{x, 0,L},
PlotRange->{-0.2,0.2}],
{t, 0, TMax, 0.1}]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
GridData=InterpolatingFunctionGrid[uSol];
ListPlot[GridData]
Dimensions[GridData]

It has taken 363 spatial dots and a total of 44 time steps of different sizes.

MMa Solver

The Mathematica solver is pretty simple at the bottom of it all.  Pick some points (equally spaced in x by default)
and semi-discretize the equations by replacing all spatial derivatives with fintie differences.  Run the ODE solver on the resulting system of ODEs.

Where it gets complicated is in knowing how many points you need for the spatial discretization and what it means about the original PDE if something goes wrong with the numerical ODE simulation.

Here is a preliminary implemtation that thinks of each of the values as a seperate function.

Clear[SemiDiscEq1,n,Δx,u,t]
SemiDiscEq1[n_,Δx_,{u_Symbol,t_Symbol}] := Module[
{Eq,x},
Eq=D[u[x,t],t]- 6 u[x,t] D[u[x,t],x] + D[u[x,t],{x,3}]==0;
Table[
Temp=Eq/.{
u[x,t]->(Subscript[u, i])[t],
D[u[x,t],t]->D[(Subscript[u, i])[t],t],
D[u[x,t],x]->(((Subscript[u, i+1])[t]-(Subscript[u, i-1])[t])/(2 Δx)),
D[u[x,t],{x,3}]->((-(Subscript[u, i-1])[t]+3 (Subscript[u, i])[t]-3 (Subscript[u, 1+i])[t]+(Subscript[u, 2+i])[t])/((Δx)^3))
};
Temp/.{
(Subscript[u, -1])->(Subscript[u, n-1]),
(Subscript[u, n])->(Subscript[u, 0]),
(Subscript[u, n+1])->(Subscript[u, 1]),
(Subscript[u, n+2])->(Subscript[u, 2])
},
{i, 0, n-1}]
]

Here is a version that thinks of the values as a vector function of time which is much easier to feed into NDSolve.

Clear[SemiDiscRHS,n,Δx,u,t,Derivative]
SemiDiscRHS[n_,Δx_,u_Symbol] := Module[
{RHS,x,t,j},
RHS= 6 u[x,t] D[u[x,t],x] - D[u[x,t],{x,3}];
Table[
Temp=RHS/.{
u[x,t]->(Subscript[u, i])[t],
D[u[x,t],x]->(((Subscript[u, i+1])[t]-(Subscript[u, i-1])[t])/(2 Δx)),     (* (Subscript[δ, 0x ])*)
(* D[u[x,t],x]->(((Subscript[u, i+1])[t]-(Subscript[u, i])[t])/(Δx)),     (* (Subscript[δ, +x ])*) *)
(* D[u[x,t],x]->(((Subscript[u, i])[t]-(Subscript[u, i-1])[t])/(Δx)),     (* (Subscript[δ, -x ])*) *)
D[u[x,t],{x,3}]->((-(Subscript[u, i-1])[t]+3 (Subscript[u, i])[t]-3 (Subscript[u, 1+i])[t]+(Subscript[u, 2+i])[t])/((Δx)^3)) (* (Subscript[δ, +x ])(Subscript[δ, 2x]) *)
(* D[u[x,t],{x,3}]->((-(Subscript[u, -2+i])[t]+3 (Subscript[u, -1+i])[t]-3 (Subscript[u, i])[t]+(Subscript[u, 1+i])[t])/((Δx)^3))       (* (Subscript[δ, -x ])(Subscript[δ, 2x]) *) *)
(* D[u[x,t],{x,3}]->((-(Subscript[u, -2+i])[t]+2 (Subscript[u, -1+i])[t]-2 (Subscript[u, 1+i])[t]+(Subscript[u, 2+i])[t])/(2((Δx)^3)))       (* (Subscript[δ, 0x ])(Subscript[δ, 2x]) *) *)
};
Temp=Temp/.{
(Subscript[u, -2])->(Subscript[u, n-2]),
(Subscript[u, -1])->(Subscript[u, n-1]),
(Subscript[u, n])->(Subscript[u, 0]),
(Subscript[u, n+1])->(Subscript[u, 1]),
(Subscript[u, n+2])->(Subscript[u, 2])
};
Simplify[Temp/.(Subscript[u, j_])[t]:>u[[j+1]]],
{i, 0, n-1}]
]

I computed the differences using the following commands e.g for everything centered.

Clear[cdiff,i,u,Δx]
cdiff[i_] := (((Subscript[u, i+1])[t]-2 (Subscript[u, i])[t]+(Subscript[u, i-1])[t])/((Δx)^2))
Simplify[((cdiff[i+1]-cdiff[i-1])/(2 Δx))]

Lets see how long it takes to work and how well it works!

Clear[L,f,TMax,n,Subscript,u]
L= 200;
IC[x_]=0.2(E^(-0.1 (((x-L/2))^2)));
TMax=12;
n=360;
Δx=N[(L/n)];
f[u_List]=SemiDiscRHS[n,Δx,u];
u0=Table[ IC[x],{x, 0.0, L-Δx, Δx}];
Timing[
uSol[t_]=u[t]/.NDSolve[{
u'[t]==f[u[t]],
u[0]==u0
},
u[t],
{t, 0, TMax}][[1]]
]

It looks like it is working OK (4 seconds for the first run) with whatever differences there are in the scheme at the moment.

Data=Table[ uSol[t],{t, 0, TMax, TMax/55}];
ListPlot3D[Data,
PlotRange->All]

Stability and Accuracy of our Scheme

Lets look at the linearized-stability for our various schemes.  Note: in this case linearized means replace 6 u (Subscript[u, x])  with 6 U (Subscript[u, x]) where U is a constant.

LinearizedRHS= 6 U D[u[x,t],x] - D[u[x,t],{x,3}]/.{
u[x,t]->(Subscript[u, j])[t],
D[u[x,t],x]->(((Subscript[u, j+1])[t]-(Subscript[u, j-1])[t])/(2 Δx)),     (* (Subscript[δ, 0x ])*)
(* D[u[x,t],x]->(((Subscript[u, j+1])[t]-(Subscript[u, j])[t])/(Δx)),     (* (Subscript[δ, +x ])*) *)
(* D[u[x,t],x]->(((Subscript[u, j])[t]-(Subscript[u, j-1])[t])/(Δx)),     (* (Subscript[δ, -x ])*) *)
D[u[x,t],{x,3}]->((-(Subscript[u, j-1])[t]+3 (Subscript[u, j])[t]-3 (Subscript[u, 1+j])[t]+(Subscript[u, 2+j])[t])/((Δx)^3)) (* (Subscript[δ, +x ])(Subscript[δ, 2x]) *)
(* D[u[x,t],{x,3}]->((-(Subscript[u, -2+j])[t]+3 (Subscript[u, -1+j])[t]-3 (Subscript[u, j])[t]+(Subscript[u, 1+j])[t])/((Δx)^3))       (* (Subscript[δ, -x ])(Subscript[δ, 2x]) *) *)
(* D[u[x,t],{x,3}]->((-(Subscript[u, -2+j])[t]+2 (Subscript[u, -1+j])[t]-2 (Subscript[u, 1+j])[t]+(Subscript[u, 2+j])[t])/(2((Δx)^3)))       (* (Subscript[δ, 0x ])(Subscript[δ, 2x]) *) *)
}

Almsot as before, substitute the k-wiggle von Neumann form of length n on the interval of length L (and the Ls do not mess anything up once we have discretized)
	(Subscript[u, j])[t,k,n]=(E^(λ t))(E^(((2π )/L) I k  (L(j/n)) ))
into our ODE.  The only difference is that we have an ODE for the amplitude of each mode and it grows like (E^(λ t)) rather than (λ^n).   Subsituting in gives
	λ (E^(λ t))(E^(2π  I k  ((j/n)) ))=(E^(λ t))α (E^(2π  I k  ((j/n)) ))
where α  is a fairly simple expression that depends on which difference we choose.  Solving gives λ=α and the mode grows if Re[α]>0 and decays if Re[λ]<0.  Here is the computation for some of the differences above.

Clear[j,n,U,k,L,Δx]
λ=Simplify[(LinearizedRHS/(E^(2π  I k  ((j/n)) )))/.(Subscript[u, j_])[t]:>(E^(2π  I k  ( (j/n)) ))]

We can play the same trick as before and think about the discrete variable ((2 k π)/n) as an angle θ and look at the behavior of all of the exponentials as trig functions.

Clear[F,θ,Δx,U]
F[θ_, Δx_, U_]=Simplify[λ/.{k->n (θ/(2 π))}]
Δx=0.1
Plot3D[ -Re[F[θ,Δx,U]], 
{θ, -π, π},
{U, -1,1}
]

Not completely clear if the thing is always stable!

Collect[λ, Δx]

It comes in two different pieces one is dominant when (U/(Δx)) is big!  the other when (1/((Δx)^3)) is big!  we actually know where each bit comes from!

We can replace the ODE y'=f(y)  with a simple ODE stepper e.g.
	({{(y^(m+1)), =, (y^m)+Δt f((y^m)), Euler}, {(y^(m+1)), =, (y^m)+Δt f((y^(m+1))), Implicit Euler}, {(y^(m+1)), =, (y^m)+((Δt)/2) (f((y^m))+f((y^(m+1)))), Heun (what we used in CN)}, {…,  ,  ,  }})
and analyze the resulting fully-discretized PDE solver with the same technique we used before.   Lets do one of them. We replace the (E^(λ t)) with (λ^m) and solve for the new λ.   Here is Euler.

Clear[j,n,U,k,L,Δx,λ,DiscreteStability]
DiscreteStability[θ_,Δt_,Δx_,U_]=Collect[
Simplify[ λ/.Solve[(λ^(m+1)) (E^(2π  I k  ((j/n)) ))==(λ^m) (E^(2π  I k  ((j/n)) )) + Δt LinearizedRHS/.(Subscript[u, j_])[t]:>(λ^m)(E^(2π  I k  ( (j/n)) )), λ][[1]]],
Δx]/.{k->n (θ/(2 π))}

Lets see what the stability plot looks like now.  Remember in the discrete case we are concerned about |λ|>1.  What deos the following mean?

{ΔxVal,UVal}={0.1, 0.1};
ContourPlot[
Abs[DiscreteStability[θ,Δt,ΔxVal,UVal]],
{θ,-π,π},
{Δt,-((ΔxVal)^3),((ΔxVal)^3)},
Contours->{0.9,1,1.1},
FrameLabel->{"θ","Δt"}
]

We can do another version here is the Trapezoid.

Clear[j,n,U,k,L,Δx,λ,DiscreteStability]
DiscreteStability[θ_,Δt_,Δx_,U_]=Collect[
Simplify[ λ/.Solve[
λ (E^(2π  I k  ((j/n)) ))== (E^(2π  I k  ((j/n)) )) + ((Δt)/2)( LinearizedRHS/.(Subscript[u, j_])[t]:>(E^(2π  I k  ( (j/n)) )) +LinearizedRHS/.(Subscript[u, j_])[t]:>λ (E^(2π  I k  ( (j/n)) ))),λ]][[1]],
Δx]/.{k->n (θ/(2 π))}

{ΔxVal,UVal}={0.1, 0.1};
ContourPlot[
Abs[DiscreteStability[θ,Δt,ΔxVal,UVal]],
{θ,-π,π},
{Δt,0,((ΔxVal)^3)},
Contours->{0.9,1,1.1},
FrameLabel->{"θ","Δt"},
RotateLabel->False
]

Break

Spectral Solver

The spectrally easy portions of the KdV
	(Subscript[u, t])+(Subscript[u, xxx])- 6 u (Subscript[u, x])=0
are the (Subscript[u, xxx]) and the (Subscript[u, x]).  Transforming we have
	(Subscript[U, t ])+𝔻3 .U - 6 U*(𝔻U)=0
where we have remembered that the transform of a product is a convolution "*" and written 𝔻 for the operater in the transform space that "performs" spatial differentiation.

Reminder about wiggle or Fourier modes

On an interval of length 2π the world is easy.

FMode[k_,n_] := Table[(E^(2 π I k (j/n))),{j,0.0,n-1}]
x=FMode[2,123];
ListPlot[{Re[x],Im[x]}]

We need to work out what to do on intervals of different lengths.

Fourier modes of derivatives.

On an interval of length 2π the world is easy.  The derivative of  y=(E^(2π I k x)) is 2π I k y in other words each frequency gets mutiplied by  2 π I k. in the derivative in a very real sense 𝔻 is a diagonal matrix with 
	(Subscript[𝔻, kk])=2 π I k
and of course (𝔻^3) is just as easily computed.  We need to work out what to do on intervals of different lengths.

Spectral Code: Version 1

(Subscript[U, t ])+𝔻3 U - 6 U*(𝔻U)=0

Clear[SpEq1, f,L, x, n,Δx,𝔻]

SpEq1=U'[t]+𝔻3.U[t] - 6 ListConvolve[U[t], 𝔻.U[t],{1,1}]==0
L=2π;
f[x_]=(E^(-(((x-L/2))^2)))
n= 123;
Δx=(L/n);
TMax=1;
𝔻 = SparseArray[ 
Band[{1,1}]->Table[ 2 π I k,{k,0.0, n-1}],
{n,n}
];
𝔻3 =  SparseArray[ 
Band[{1,1}]->Table[ (((2 π I k))^3),{k,0.0, n-1}],
{n,n}
];
u0=Table[ f[x], {x, 0.0, L- Δx, Δx}];
U0=Fourier[u0];
USol[t_]=U[t]/.NDSolve[{
U'[t]+𝔻3.U[t] - 6 ListConvolve[U[t], 𝔻.U[t],{1,1}]==0,
U[0]==U0
},
U[t],
{t, 0, TMax},
Compiled->False,
MaxSteps->(10^5),
Method->"BDF"][[1]]

Manipulate[ 
ListPlot[ Chop[InverseFourier[ USol[t]]]],
{t, 0, 0.4,0.01}]

Misc

Do we need all the Fourier modes.

One option is to keep all the wiggle modes.  The other is to only keep some number of the lower modes" this brings up the question of what to do with the ones you are going to neglect.