ODE Test Problems

We need a variety of test problems in different dimensions etc.  The generic format is 
	y'=f[y]
and all we really need is a formula for f:(ℝ^n)->(ℝ^n) and we have a test problem.  We want to make a range of problems in different sizes and types.

Small named problems.

We have a collection of known problems.  SHM, Orbit, van der Pol, Brusselator, Robertson-Stiff, Duffing, Swing-Spring, Mandelbrot, Fitzhugh-Nagumo, Hodgkin-Huxley etc.

van der Pol

The second order form.

μ=3;
TMax=15;
xSol[t_]=x[t]/.NDSolve[{ 
x''[t]-μ (1 - ((x[t])^2)) x'[t]+x[t]==0,
x[0]==1,
x'[0]==2
},
x[t],
{t,0, TMax}][[1]];
{Plot[ xSol[t],{t, 0,TMax}] ,ParametricPlot[{xSol[t],xSol'[t]},{t, 0, TMax},
PlotRange->All]}

A slighlt different scalin is going on in the second order form on wiki.

μ=3;
TMax=15;
fVDP[{x1_,x2_}]={
μ ( x1- (1/3) (x1^3)- x2),
(1/(μ ))x1}
xSol[t_]=x[t]/.NDSolve[{ 
x'[t]==fVDP[x[t]],
x[0]=={1,2}
},
x[t],
{t,0, TMax}][[1]];
{Plot[ xSol[t][[1]],{t, 0,TMax}] ,ParametricPlot[xSol[t],{t, 0, TMax},
PlotRange->All]}

Extending small problems

We can embed small problems in to linear problems with controlled eigenvalue distributions. One way to do this is to have a clean system and then roatet the coordinates.  To do this I need to build a rotated coordinate frame. 
	x'=f[x
mulitply by a fixed rotation/reflection Q
	Q x'=Q f[x]
then combine the derivative with the Q 
	  (Q x)'=Q f[(Q^T)Q x]
and define a new variable y= Q x to get
	y'=Q f[(Q^T)y]
which will look all messed up but will actually just be a rotated version of the original.

RandomQ[n_]:=Module[
{Q},
Q=QRDecomposition[RandomReal[{-1,1},{n,n}]][[1]]
]

fVDPBig[ {x1_,x2_,x3_,x4_,x5_}]={
μ ( x1- (1/3) (x1^3)- x2),
(1/(μ ))x1,
- x3,
-2 x4,
-0.1 x5
};
Q=RandomQ[5];
fVDPBigMess[y_]:=Q.fVDPBig[Transpose[Q].y];

μ=3;
TMax=100;
xSol[t_]=x[t]/.NDSolve[{
x'[t]==fVDPBig[x[t]],
x[0]=={1,2,3,4,5}
},
x[t],
{t,0, TMax}][[1]]
ySol[t_]=y[t]/.NDSolve[{
y'[t]==fVDPBigMess[y[t]],
y[0]==Q.{1,2,3,4,5}
},
y[t],
{t,0, TMax}][[1]]
zSol[t_]:=Transpose[Q].ySol[t]
{ParametricPlot[{xSol[t][[{1,2}]],ySol[t][[{1,2}]]},
{t, 0, TMax}],
ParametricPlot[{xSol[t][[{1,2}]],zSol[t][[{1,2}]]},
{t, 0, TMax}]}

Method of Lines

A PDE generates a large system of ODEs by replacing all the spatial derivatives with finite difference aproximations.  These tend to have specific eigenvalue distributions for their linear parts.

As an example w will do the heat equation
	(Subscript[u, t])=α (Subscript[u, xx])+f[x]
on the interval of length π with periodic boundary conditions
	u[0,t]=u[π,t]
and a generic intial condition
	u[x,0]=g[x].
We pick a number of points "nx" discretize in x  i.e. (Subscript[x, j])=j Δx where Δx=π/n and approximate the second derivative with the difference quotient.
	(Subscript[u, xx])((Subscript[x, j]),t)≈((u[(Subscript[x, j+1]),t]-2u[(Subscript[x, j]),t]+u[(Subscript[x, j-1]),t])/((Δx)^2))
to get the system of ODEs
	(Subscript[u, j])'=a ((Subscript[u, j+1])-2 (Subscript[u, j])+(Subscript[u, j-1]))+(Subscript[f, j])
where (Subscript[u, j])[t]≈u[(Subscript[x, j]),t] and we definea=α/((Δx)^2) along with (Subscript[f, j])[t]=f[(Subscript[x, j]),t].  We implement the periodic boundary conditions by identifying (Subscript[u, 0])[t]=(Subscript[u, nx])[t]and (Subscript[u, 1])[t]=(Subscript[u, nx+1])[t] and (Subscript[u, -1])[t]=(Subscript[u, nx-1])[t]. We get the following system of ODEs.

TMax=100;
nx=12;
α=0.1;
f[x_,t_]=t (E^(-0.1t ))Sin[x];
g[x_]= Sin[3x ];
x[j_]=0.0+ j π/nx; 
Eqs=Join[
Table[
D[u[j][t],t]== α(u[j+1][t]-2u[j][t]+u[j-1][t]) + f[x[j],t],
{j,1, nx}],
Table[
u[j][0]==g[x[j]],
{j,1, nx}]
]/.{
u[0]->u[nx],
u[nx+1]->u[1]
};

uSol[t_]=Table[u[j][t],{j,1,nx}]/.NDSolve[
Eqs,
Table[u[j],{j,1,nx}],
{t, 0, TMax}
][[1]];

Plot[uSol[t],{t, 0, TMax},
AxesLabel->{"t","u[j]s"}]
ListPlot3D[
Table[uSol[t],{t, 0, TMax, TMax/40.0}],
AxesLabel->{"x","t","u"},
DataRange->{{0,π},{0, TMax}}
]

Spectral and Other Expansions

A PDE can generate a large system of ODEs via "spectral" expansions.  We will do the heat equation example again.

As an example w will do the heat equation
	(Subscript[u, t])=α (Subscript[u, xx])+f[x]
on the interval of length π with Dircihlet boundary conditions
	u[0,t]=u[π,t]=0
and a generic intial condition
	u[x,0]=g[x].
We pick a number of functions say
	(Subscript[ϕ, i])[x]  for i=1…n
which satisfy the boundary conditions and form a finite sum with time varying coefficients
	(Subscript[c, 1])[t] (Subscript[ϕ, 1])[x] + (Subscript[c, 2])[t] (Subscript[ϕ, 2])[x]+…+(Subscript[c, n])[t](Subscript[ϕ, n])[x]
and substitute this into the PDE to get
	(∑_(i=1)^m)c[m]'[t](Subscript[ϕ, i])[x]=(∑_(i=1)^m)-a c[m][t](Subscript[ϕ, i])''[x]+f[x]
and substitue into the initial condition to get
	(∑_(i=1)^m)c[m][0](Subscript[ϕ, i])[x]=g[x]

We need to get a system of ODES the c[i].  One standard trick for the heat equation is to use (Subscript[ϕ, i])[x]=Sin[ix]Sin[i x] and separate the equations by multiplying the equation by a specific (Subscript[ϕ, j])[x] and integrating fromx=0->π;  the result is a Fourier spectral expansion.  A second standard trick is to use tent functions for the (Subscript[ϕ, i])[x]  and simply evaluate the equation at the peaks of each tent;  this is a collocation technique.  You can use other shape functions and more elaborate weighted collocation or quadrature schemes.

Spectral Expansion

With (Subscript[ϕ, i])[x]=Sin[i x] and the intgration trick you get the diagonal system 
	(Subscript[c, i])'[t]=-a (i^2)(Subscript[c, i])[t]+(Subscript[f, i])[t]
	(Subscript[c, i])[0]=(Subscript[g, i])	
where 
	(Subscript[f, i])[t]=(∫_0^π)f[x,t]Sin[i x]dx
and
	(Subscript[g, i])=(∫_0^π)g[x]Sin[i x]dx.

TMax=100;
nx=12;
a=0.1;
f[x_,t_]=t (E^(-0.1t ))Sin[x];
f[1][t]=t (E^(-0.1t ));
f[2][t]=0;
f[3][t]=0;
g[x_]= Sin[3x ];
g[1]=0;
g[2]=0;
g[3]=1;
Eqs=Join[
Table[
D[c[j][t],t]== -a (j^2)c[j][t]+f[j][t],
{j,1, 3}],
Table[
c[j][0]==g[j],
{j,1, 3}]
];
cSol[t_]=Table[c[j][t],{j,1,3}]/.NDSolve[
Eqs,
Table[c[j],{j,1,3}],
{t, 0, TMax}
][[1]];

Plot[cSol[t],{t, 0, TMax},
AxesLabel->{"t","u[j]s"}]

Plot3D[
Sum[ cSol[t][[j]] Sin[j x],{j,1,3}],
{x,0,π},
{t,0, TMax}]

Tent Expansion

Lets pick a bunhc of equally spaced points.
	Δx=π/nx
	(Subscript[p, i])=i Δx  for i=0,…,nx
and define tent functions centered on each interior point.

nx=6;
Δx=π/nx;
p[i_]= i Δx;
ϕ[i_][x_]= Piecewise[{
{0,x ≤ (i-1) Δx},
{1-i+(x/(Δx)),(i-1) Δx<x≤i Δx},
{1+i-(x/(Δx)),i  Δx<x≤(i+1) Δx},
{0,x> (i+1)Δx}
}
];
Plot[ Table[ϕ[i][x],{i,1,nx-1}],{x,0,π}]
f[x_]= Sin[x];
fApp[x_]=Sum[f[i Δx] ϕ[i][x],{i,1,nx-1}];
Plot[{f[x],fApp[x]},{x,0,π}]

Large scale physical problems e.g. Saturn rings and moons.

Include a large number of interacting physical objects. e,g  gravitationally build the model automatically.