4620 Homework 3

Utilities

Ap1

Ap1[n_]:=SparseArray[{
Band[{1,1}]->-2.0,
Band[{1,2}]->  1.0,
Band[{2,1}]->  1.0,
{1,n}->1.0,
{n,1}->1.0
},
{n,n}]
MatrixForm[Ap1[6]]

Ap2

Ap2[n_]:=SparseArray[{
Band[{1,1}]->-2.0,
Band[{1,2}]->  1.0,
Band[{2,1}]->  1.0,
{1,n-1}->1.0,
{n,2}->1.0
},
{n,n}]
MatrixForm[Ap2[6]]

Md

Md[n_,nu_]:=SparseArray[{
Band[{2,2},{n-1,n}]->1-2.0nu,
Band[{2,3},{n-1,n}]->  nu,
Band[{2,1},{n-1,n}]->  nu
},
{n,n}]
MatrixForm[Md[6,nu]]

Solutions

{a,b}={0.0,2π};
uExact[x_,t_]=2+4 (E^(-(3^2)t))Cos[3 x] + 2 (E^(-t))Sin[x]; 
u0[x_]=uExact[x,0];
ContourPlot[ uExact[x,t],{x,a,b},{t,0, 0.2}]
Plot3D[ uExact[x,t],{x,a,b},{t,0, 0.2}]

Here is our numerical solution

nx=8;
{nt,Δt}={50,0.01};
Δx=(b-a)/(nx+1);
ν=Δt/((Δx)^2)
Mp=IdentityMatrix[nx+2]+ν Ap2[nx+2];
USol=ConstantArray[0,{nt,nx+2}];     (* storage "a" *)
USol[[1]] = Table[u0[x],{x,a,b,Δx}]; (* initial condition "b"*) 
Do[                                                                            (* Iteration *)
USol[[n+1]]=Mp.USol[[n]],
{n,1,nt-1}]
UExact=Table[uExact[x,i Δt],{i,0,nt-1},{x,a,b,Δx}];
((Norm[USol-UExact])/(Norm[UExact]))

Manipulate[
Show[
Plot[uExact[x, i Δt],{x,a,b},
PlotRange->{-5.5,5.5}],
ListPlot[USol[[i+1]],
DataRange->{a,b}]
],
{i,0, nt,1}]

{a,b}={0.0,2π};
uExact[x_,t_]=4 (E^(-25t))Sin[5 x] + 2 (E^(-9 t))Sin[3x]; 
u0[x_]=uExact[x,0];
ContourPlot[ uExact[x,t],{x,a,b},{t,0, 0.2}]
Plot3D[ uExact[x,t],{x,a,b},{t,0, 0.2}]

nx=32;
{nt,Δt}={50,0.01};
Δx=(b-a)/(nx+1);
ν=Δt/((Δx)^2)
MdMat=Md[nx+2,ν];
USol=ConstantArray[0,{nt,nx+2}];     (* storage "a" *)
USol[[1]] = Table[u0[x],{x,a,b,Δx}]; (* initial condition "b"*) 
Do[                                                                            (* Iteration *)
USol[[n+1]]=MdMat.USol[[n]],
{n,1,nt-1}]

ListPlot3D[USol,
PlotRange->All]

UExact=Table[uExact[x,i Δt],{i,0,nt-1},{x,a,b,Δx}];
((Norm[USol-UExact])/(Norm[UExact]))

Manipulate[
Show[
Plot[uExact[x, i Δt],{x,a,b},
PlotRange->{-5.5,5.5}],
ListPlot[USol[[i+1]],
DataRange->{a,b}]
],
{i,0, nt,1}]

Not saviong all the data!

nx=50;
{nt,Δt}={50,0.01};
Δx=(b-a)/(nx+1);
ν=Δt/((Δx)^2)
MdMat=Md[nx+2,ν];
u = Table[u0[x],{x,a,b,Δx}]; (* initial condition "b"*) 
Do[                                                                            (* Iteration *)
u=MdMat.u,
{nt-1}]
Show[
ListPlot[u,
DataRange->{a,b}],
Plot[uExact[x,nt*Δt],{x,a,b}]
]