Example 1 Wave Equation Dirichlet Conditions 

Assignment of system parameters  

>

restart:

 

>	L:=Pi; c:=1;

 

>	with(plots):

 

>

 

These notes are concerned with solving the heat equation using Fourier's method 

First we consider Dirichlet boundary conditions:   

>

 

In this case we obtain the following eigenvalues and orthonormal eigenfunctions.  

>	mu[n]:=(n*Pi/L); lambda[n]:=-mu[n]^2;

 

>

 

for n = 1, 2, 3, ...   Orthonormal eigenfunctions 

>	X[n](x):=sqrt(2/L)*sin(n*Pi/L*x);X[m](x):=subs(n=m,X[n](x)):

 

>

 

Statement of orthonormality  

>	Int(X[n](x)*X[m](x),x=0...L)=delta(n,m);

 

Time-dependent solution  

>	T[n](t):=a(n)*cos(c*mu[n]*t)+b(n)*sin(c*mu[n]*t);

 

>

 

>	u[n](x,t):=T[n](t)*X[n](x);

 

>

 

>

 

Eigenfunction expansion  

>	u(x,t):=Sum(u[n](x,t),n=1..infinity);

 

Evaluation of Fourier coefficients for the specific initial condition  

>	f(x):=x*(Pi-x);

 

>	a(n):=Int(f(x)*X[n](x),x=0..L);a(n):=expand(value(%)):

 

>

 

>

 

>	a(n):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},a(n)));

 

>	g(x):=x^2*(Pi-x);

 

>	b(n):=1/(mu[n]*c)*Int(g(x)*X[n](x),x=0..L);b(n):=expand(value(%)):

 

>	b(n):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},b(n)));

 

>

 

>

 

>

 

>

 

Check Fourier coefficients  

>	plot({f(x),sum(a(n)*X[n](x),n=1..3)},x=0..L);

 

>

 

>

 

>	plot({g(x),sum(mu[n]*c*b(n)*X[n](x),n=1..5)},x=0..L);

 

>

 

>

 

Our simple solutions 

>	u[n](x,t):=eval(T[n](t)*X[n](x));

 

Series solution 

>	u(x,t):=Sum(u[n](x,t),n=1..infinity);

 

First few terms of sum 

>	u(x,t):=sum(u[n](x,t),n=1..3);

 

Animation 

>	animate(plot,[u(x,t),x=0..L],t=0..15,color=red,thickness=3);

 

>

 

Animation Sequence 

>	u(x,0):=subs(t=0,u(x,t)):u(x,1):=subs(t=1,u(x,t)):

 

>	u(x,2):=subs(t=2,u(x,t)):u(x,3):=subs(t=3,u(x,t)):

 

>	u(x,4):=subs(t=4,u(x,t)):u(x,5):=subs(t=5,u(x,t)):

 

>	plot({u(x,0),u(x,1),u(x,2),u(x,3),u(x,4),u(x,5)},x=0..L,thickness=3);

 

>	plot3d(u(x,t),x=0..L,t=0..5,axes=boxed);

 

>

 

Example 2 Wave Equation Neumann Conditions 

Assignment of system parameters  

>

restart:

 

>	L:=Pi; c:=1;

 

>	with(plots):

 

In this section  we consider Neumann  boundary conditions:   

>

 

>

 

In this case the eigenvalues and orthonormal eigenfunctions are given by 

>	mu[n]:=(n*Pi/L); lambda[n]:=-mu[n]^2; lambda[0]:=0;

 

>

 

for n = 1, 2, 3, ...   Orthonormal eigenfunctions 

>	X[n](x):=sqrt(2/L)*cos(n*Pi/L*x);X[m](x):=subs(n=m,X[n](x)):

 

>	X[0](x):=1/sqrt(L);

 

Statement of orthonormality  

>	Int(X[n](x)*X[m](x),x=0...L)=delta(n,m);

 

Time-dependent solution  

>	T[n](t):=a(n)*cos(c*n*t)+b(n)*sin(c*n*t);

 

>

 

>	u[n](x,t):=T[n](t)*X[n](x);

 

>	u[0](x,t):=(a(0)+b(0)*t)/2*X[0](x);

 

>

 

Eigenfunction expansion  

>	u(x,t):=u[0](x,t)+Sum(u[n](x,t),n=1..infinity);

 

Evaluation of Fourier coefficients for the specific initial condition  

>	f(x):=x*(Pi-x);

 

>

 

>	a(0):=Int(f(x)*X[0](x) ,x=0..L);a(0):=expand(value(%)):

 

>	a(0):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},a(0)));

 

>	a(n):=Int(f(x)*X[n](x),x=0..L);a(n):=expand(value(%)):

 

>	a(n):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},a(n)));

 

>

 

Check Fourier coefficients  

>	plot({f(x), a(0)*X[0](x)+sum(a(n)*X[n](x),n=1..7)},x=0..L);

 

>

 

>	g(x) := 0;

 

>

 

>	b(0):=Int(g(x)*X[0](x),x=0..L);b(0):=expand(value(%)):

 

>	b(0):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},b(0)));

 

>

 

>	b(n):=Int(g(x)*X[n](x),x=0..L);b(n):=expand(value(%)):

 

>	b(n):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},b(n)));

 

>

 

>

 

Check Fourier coefficients  

>	plot({g(x), b(0)*X[0](x)+sum(b(n)*X[n](x),n=1..7)},x=0..L);

 

>

 

Our simple solutions 

>	u[n](x,t):=eval(T[n](t)*X[n](x));

 

>	u[0](x,t):=(a(0)+b(0)*t)/2*X[0](x);

 

Series solution 

>	u(x,t):=u[0](x,t)+Sum(u[n](x,t),n=0..infinity);

 

First few terms of sum 

>	u(x,t):=u[0](x,t)+sum(u[n](x,t),n=1..3);

 

Animation 

>	animate(plot,[u(x,t),x=0..L],t=0..15,color=red,thickness=3);

 

Animation Sequence 

>	u(x,0):=subs(t=0,u(x,t)):u(x,1):=subs(t=1,u(x,t)):

 

>	u(x,2):=subs(t=2,u(x,t)):u(x,3):=subs(t=3,u(x,t)):

 

>	u(x,4):=subs(t=4,u(x,t)):u(x,5):=subs(t=5,u(x,t)):

 

>	plot({u(x,0),u(x,1),u(x,2),u(x,3),u(x,4),u(x,5)},x=0..L,thickness=3);

 

>	plot3d(u(x,t),x=0..L,t=0..5,axes=boxed);

 

>

 

Example 3 Wave Equation Dirichlet-Neumann Conditions 

Assignment of system parameters  

>

restart:

 

>	L:=1; c:=1;

 

>	with(plots):

 

>

 

In this section  we consider Dirchlet and Neumann  boundary conditions:   

>

 

>

 

In this case the eigenvalues and orthonormal eigenfunctions are given by  

>	mu[n]:=((2*n-1)*Pi/(2*L)); lambda[n]:=-mu[n]^2;

 

>

 

for n = 1, 2, 3, ...   Orthonormal eigenfunctions 

>	X[n](x):=sqrt(2/L)*sin(mu[n]*x);X[m](x):=subs(n=m,X[n](x)):

 

>

 

Statement of orthonormality  

>	Int(X[n](x)*X[m](x),x=0...L)=delta(n,m);

 

>	radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},int(X[n](x)^2,x=0..L)));

 

Time-dependent solution  

>	T[n](t):=a(n)*cos(c*n*t)+b(n)*sin(c*n*t);

 

>

 

>	u[n](x,t):=T[n](t)*X[n](x);

 

>

 

>

 

Eigenfunction expansion  

>	u(x,t):=Sum(u[n](x,t),n=1..infinity);

 

Evaluation of Fourier coefficients for the specific initial condition  

>	f(x):=x*sin(Pi*x);

 

>	a(n):=Int(f(x)*X[n](x),x=0..L);a(n):=expand(value(%)):

 

>	a(n):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},a(n)));

 

>

 

Check Fourier coefficients  

>	plot({f(x),sum(a(n)*X[n](x),n=1..5)},x=0..L);

 

>

 

>

g(x):=0;

 

>	b(n):=1/(n*c)*Int(g(x)*X[n](x),x=0..L);b(n):=expand(value(%)):

 

>	b(n):=radsimp(subs({sin(n*Pi)=0,cos(n*Pi)=(-1)^n},b(n)));

 

>

 

Check Fourier coefficients  

>	plot({g(x),sum(b(n)*X[n](x),n=1..5)},x=0..L);

 

>

 

>

 

>

 

>

 

Our simple solutions 

>	u[n](x,t):=eval(T[n](t)*X[n](x));

 

Series solution 

>	u(x,t):=Sum(u[n](x,t),n=1..infinity);

 

First few terms of sum 

>	u(x,t):=sum(u[n](x,t),n=1..3);

 

Animation 

>	animate(plot,[u(x,t),x=0..L],t=0..15,color=red,thickness=3);

 

Animation Sequence 

>	u(x,0):=subs(t=0,u(x,t)):u(x,1):=subs(t=1,u(x,t)):

 

>	u(x,2):=subs(t=2,u(x,t)):u(x,3):=subs(t=3,u(x,t)):

 

>	u(x,4):=subs(t=4,u(x,t)):u(x,5):=subs(t=5,u(x,t)):

 

>	plot({u(x,0),u(x,1),u(x,2),u(x,3),u(x,4),u(x,5)},x=0..L,thickness=3);

 

>	plot3d(u(x,t),x=0..L,t=0..5,axes=boxed);

 

>