Fall 2012

Code Complexity

Students find the code very hard and once it gets too complicated for them they shut down.
Students do not get swamped with the Sparse Array structure.  They like band.
They have no problem with periodic boundary conditions. 
They have a terrible time (and I do as well) with the n subsidivisions and n+1 points.
They do get the implicit/explicit dichotomy.  With a sparse linear solve it becomes less clear that one is harder than the other.
They have issues with the Thomas algorithm for periodic stuff.  Not sure how to answer this since
they have real problems with any BC except periodic.
Dirichlet has problems switiching explicit/implicit

Analytical Complexity

Students find LTE OK.  But it needs to be very clear and explicit what you do.
Students find the remainder terms hard particularly the floating point form.
Students like when it goes back to calc 1
Students are fine with a library of differences but it should be short.
Students are fine with errors for these differences but again the library should be short.
Software support for series and trig stuff is good.

Fourier Complexity

Students have a very spotty exposure to Fourier Series.  
The transition to discrete is problematic.
The highest resolvvable mode is tricky for them to get.
The mismatch between (0,1), (-π,π), (0,2π) etc. is confusing for them.
There are too many things in the Morton and Mayer discussion.
The ±∞ in the dicrete sum is tricky for students.  Bth because of the ∞ and the ±. 
The highest resolvvable mode is very tricky for them to get. 
Good entree to spectral solvers.

Computational Complexity

The sparse Array and solver are a blessing.
The sparse linear solve is a cool way to motivate iterative solvers and/or the Thomas Algorithm.
Also a good proframming exercise.

1D versus multi dimensional problems

The book-keeping for multi-D is very hard for students. 
It is even harder to grade. 
1D first.

Sample Code

(Subscript[u, t])=k (Subscript[u, xx])+f[x,t]  with homogeneous Dirichlet conditions.  We will use test problem below to check the accuracy of our code.

Clear[a,b,u0,x]
{a,b}={0.0,1.0};
nx=48;
Δx = (((b-a))/nx);
u0[x_]= Sin[π x] + 3 Sin[2 π x]-4 Sin[4 π x];
U0=Table[u0[x],{x,a,b,Δx}];
Show[
Plot[u0[x],{x,a,b}],
ListPlot[U0,
DataRange->{a,b}],
PlotLabel->"IC"
]

TMax=2;
f[x_,t_]= Cos[t]Sin[π x] - 3Sin[t] Sin[ 3 π x] ;
ContourPlot[f[x,t],{x,a,b},{t,0,TMax}]

The analytic solution (with f=0) from separation of variables is

uh[x_,t_]=  (E^(- (((π))^2)t))Sin[π x] + 3 (E^(-(((2π))^2)t)) Sin[2 π x]-4 (E^(-(((4 π))^2)))Sin[4 π x];
ContourPlot[uh[x,t],{x,a,b},{t,0,0.1},
PlotRange->All]

The analytic solution (with the given f) from a Fourier expansion is

{b1[t]/.DSolve[ {
b1'[t] +(π^2) b1[t]==Cos[t],
b1[0]==1},
b1[t],
t][[1]],
b2[t]/.DSolve[ {
b2'[t] +(((2π))^2) b2[t]==0,
b2[0]==3},
b2[t],
t][[1]],
b3[t]/.DSolve[ {
b3[t] +(((3π))^2) b3[t]==-3Sin[t],
b3[0]==0},
b3[t],
t][[1]],
b4[t]/.DSolve[ {
b4'[t] +(((4π))^2) b4[t]==0,
b4[0]==-4},
b4[t],
t][[1]]
}//Simplify

ue[x_,t_]={(((E^(-(π^2) t)) (1-(π^2)+(π^4)+(E^((π^2) t)) (π^2) Cos[t]+(E^((π^2) t)) Sin[t]))/(1+(π^4))),3 (E^(-4 (π^2) t)),-((3 Sin[t])/(1+9 (π^2))),-4 (E^(-16 (π^2) t))}.{Sin[π x], Sin[2 π x], Sin[3 π x], Sin[4 π x]}

ContourPlot[ue[x,t],{x,a,b},{t,0,TMax},
PlotRange->All,
PlotPoints->{40,300}
]

Wandering Sub-Domains

I need to make an image of the rectangular 2D domain.

{nx,ny}={128,64}
dim=16;
Manipulate[
Show[
Graphics[{
{LightPink,Rectangle[{0,0},{nx,ny}]},
{Blue,Opacity[0.5],Rectangle[p-Floor[dim/2],p+Floor[dim/2]]}
}
],
Frame->True,
PlotRange->All,
AspectRatio->Automatic,
PlotLabel->"Single Computational Sub Domain"
],
{{p,Floor[{nx,ny}/2]},Locator}
]

Multiple subdomains

{nx,ny}={128,64}
dim=16;
Manipulate[
Show[
Graphics[{
{LightPink,Rectangle[{0,0},{nx,ny}]},
{Blue,Opacity[0.5],Rectangle[p1-Floor[dim/2],p1+Floor[dim/2]]},
{Green,Opacity[0.5],Rectangle[p2-Floor[dim/2],p2+Floor[dim/2]]},
{Red,Opacity[0.5],Rectangle[p3-Floor[dim/2],p3+Floor[dim/2]]},
{Cyan,Opacity[0.5],Rectangle[p4-Floor[dim/2],p4+Floor[dim/2]]}
}
],
Frame->True,
PlotRange->All,
AspectRatio->Automatic,
PlotLabel->"Muliple Computational Sub Domain"
],
{{p1,Floor[{nx/16,7ny/8}]},Locator},
{{p2,Floor[{nx/16,5ny/8}]},Locator},
{{p3,Floor[{nx/16,3ny/8}]},Locator},
{{p4,Floor[{nx/16,ny/8}]},Locator}
]

You can overlap domains if you want!  You can move each of them seperately if you want.  You can focus on one particular region fi you would like.   In the CUDA world version each moves independently.  You can ensure that writes happen in each step.   You can move things in a number of different ways.   Back to a single domain to illustrate the algorithm for a particular subdomain motion.

{nx,ny}={128,64}
dim=16;
Manipulate[
Show[
p={n,Floor[ny/2]};
Graphics[{
{LightPink,Rectangle[{0,0},{nx,ny}]},
{Blue,Opacity[0.5],Rectangle[p-Floor[dim/2],p+Floor[dim/2]]},
p=p+{3,0};
{Red,Opacity[0.5],Rectangle[p-Floor[dim/2],p+Floor[dim/2]]}
}
],
Frame->True,
PlotRange->All,
AspectRatio->Automatic,
PlotLabel->"Single Computational Sub Domain"
],
{n, dim/2, nx-1-dim/2,1}
]

Read in Red, compute in purple; write in blue. The border cells can only be partially relaxed.

Relaxation Iteration

The relaxation iteration is basically a potentially asynchronus  SOR iteration with the known optimal relaxation parameter from the PDE. The Crank-Nicolson scheme is 
	(((Subscript[𝒰, i,j]^(n+1))-(Subscript[𝒰, i,j]^n))/(Δt))=(k/(2(h^2)))( (({{ , (Subscript[𝒰, i,+1j]^(n+1)),  }, {(Subscript[𝒰, i,j-1]^(n+1)), -4 (Subscript[𝒰, i,j]^(n+1)), (Subscript[𝒰, i,j+1]^(n+1))}, { , (Subscript[𝒰, i-1,j]^(n+1)),  }}))+ (({{ , (Subscript[𝒰, i,+1j]^n),  }, {(Subscript[𝒰, i,j-1]^n), -4 (Subscript[𝒰, i,j]^n), (Subscript[𝒰, i,j+1]^n)}, { , (Subscript[𝒰, i-1,j]^n),  }})))+((((Subscript[ℱ, i,j]^(n+1))+(Subscript[ℱ, i,j]^n)))/2)
or equvalently
	(Subscript[𝒰, i,j]^(n+1))-ν  (({{ , (Subscript[𝒰, i,+1j]^(n+1)),  }, {(Subscript[𝒰, i,j-1]^(n+1)), -4 (Subscript[𝒰, i,j]^(n+1)), (Subscript[𝒰, i,j+1]^(n+1))}, { , (Subscript[𝒰, i-1,j]^(n+1)),  }}))=(Subscript[𝒰, i,j]^n)+ν (({{ , (Subscript[𝒰, i,+1j]^n),  }, {(Subscript[𝒰, i,j-1]^n), -4 (Subscript[𝒰, i,j]^n), (Subscript[𝒰, i,j+1]^n)}, { , (Subscript[𝒰, i-1,j]^n),  }}))+((Δt)/2)((Subscript[ℱ, i,j]^(n+1))+(Subscript[ℱ, i,j]^n))
which can be considered as the standard linear system for the vector  (𝕌^n)⟷(𝒰^n) where (Subscript[𝒰, i,j])⟷(Subscript[𝕌, j+i*m]) if we do zero offst vectors and matrices ainthe style of a C program.

{nx,ny}={128,64}
dim=16;
Manipulate[
Show[
Graphics[{
{LightPink,Rectangle[{0,0},{nx,ny}]},
{Blue,Opacity[0.5],Rectangle[p-Floor[dim/2],p+Floor[dim/2]]}
}
],
Frame->True,
PlotRange->All,
AspectRatio->Automatic,
PlotLabel->"Single Computational Sub Domain"
],
{{p,Floor[{nx,ny}/2]},Locator}
]

?Manipulate

Fall 2013

Computing issues dominated.
Deall ii is too complicated for the students.
Look at FENiCS.