m4354_heat_examples.mw

 

Some Examples Using Maple to Solve the Heat Equation 

 

Example 1 Heat Equation Dirichlet Conditions 

 

Assignment of system parameters  

 

> restart:
 

> L:=Pi; k:=1;
 

 

Pi
1
 

> with(plots):
 

>
 

Eigenvalues and orthonormal eigenfunctions.  

 

 

 

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

 

n
`+`(`-`(`*`(`^`(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)):
 

 

`/`(`*`(`^`(2, `/`(1, 2)), `*`(sin(`*`(n, `*`(x))))), `*`(`^`(Pi, `/`(1, 2))))
 

>
 

Statement of orthonormality  

 

 

 

 

 

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

 

Int(`+`(`/`(`*`(2, `*`(sin(`*`(n, `*`(x))), `*`(sin(`*`(m, `*`(x)))))), `*`(Pi))), x = 0 .. Pi) = delta(n, m)
 

Time-dependent solution  

 

 

> T[n](t):=C(n)*exp(k*lambda[n]*t);
 

`*`(C(n), `*`(exp(`+`(`-`(`*`(`^`(n, 2), `*`(t)))))))
 

>
 

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

`/`(`*`(C(n), `*`(exp(`+`(`-`(`*`(`^`(n, 2), `*`(t))))), `*`(`^`(2, `/`(1, 2)), `*`(sin(`*`(n, `*`(x))))))), `*`(`^`(Pi, `/`(1, 2))))
 

>
 

>
 

Eigenfunction expansion  

 

 

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

Sum(`/`(`*`(C(n), `*`(exp(`+`(`-`(`*`(`^`(n, 2), `*`(t))))), `*`(`^`(2, `/`(1, 2)), `*`(sin(`*`(n, `*`(x))))))), `*`(`^`(Pi, `/`(1, 2)))), n = 1 .. infinity)
 

 

Evaluation of Fourier coefficients for the specific initial condition  

 

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

`*`(x, `*`(`+`(Pi, `-`(x))))
 

 

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

 

Int(`/`(`*`(x, `*`(`+`(Pi, `-`(x)), `*`(`^`(2, `/`(1, 2)), `*`(sin(`*`(n, `*`(x))))))), `*`(`^`(Pi, `/`(1, 2)))), x = 0 .. Pi)
 

 

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

 

`+`(`-`(`/`(`*`(2, `*`(`^`(2, `/`(1, 2)), `*`(`+`(`-`(1), `^`(-1, n))))), `*`(`^`(Pi, `/`(1, 2)), `*`(`^`(n, 3))))))
 

>
 

Check Fourier coefficients  

 

 

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

Plot_2d
 

>
 

 

Our simple solutions 

 

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

`+`(`-`(`/`(`*`(4, `*`(`+`(`-`(1), `^`(-1, n)), `*`(exp(`+`(`-`(`*`(`^`(n, 2), `*`(t))))), `*`(sin(`*`(n, `*`(x))))))), `*`(Pi, `*`(`^`(n, 3))))))
 

 

Series solution 

 

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

 

Sum(`+`(`-`(`/`(`*`(4, `*`(`+`(`-`(1), `^`(-1, n)), `*`(exp(`+`(`-`(`*`(`^`(n, 2), `*`(t))))), `*`(sin(`*`(n, `*`(x))))))), `*`(Pi, `*`(`^`(n, 3)))))), n = 1 .. infinity)
 

 

First few terms of sum 

 

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

 

`+`(`/`(`*`(8, `*`(exp(`+`(`-`(t))), `*`(sin(x)))), `*`(Pi)), `/`(`*`(`/`(8, 27), `*`(exp(`+`(`-`(`*`(9, `*`(t))))), `*`(sin(`+`(`*`(3, `*`(x))))))), `*`(Pi)))
 

Animation 

 

> animate(plot,[u(x,t),x=0..L],t=0..10,color=red,thickness=3, numpoints=200, frames=30);
 

 

 

 

Plot_2d
 

The preceding animation command illustrates the spatial-time dependent solution for u(x, t). The animation sequence  shows snapshots at times t = 0, 1, 2, 3, 4, and 5. 

 

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);
 

Plot_2d
 

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

Plot_2d
 

>
 

 

Example 2 Heat Equation Neumann Conditions 

 

Assignment of system parameters  

 

> restart:
 

> L:=2; k:=1/10;
 

 

2
`/`(1, 10)
 

> with(plots):
 

>
 

Eigenvalues and orthonormal eigenfunctions.  

 

 

 

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

 

`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi))))
`+`(`-`(`*`(`/`(1, 4), `*`(`^`(n, 2), `*`(`^`(Pi, 2))))))
 

>
 

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)):
 

 

cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x))))))
 

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

`+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)))))
 

Statement of orthonormality  

 

 

 

 

 

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

 

Int(`*`(cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x)))))), `*`(cos(`+`(`*`(`/`(1, 2), `*`(m, `*`(Pi, `*`(x)))))))), x = 0 .. 2) = delta(n, m)
 

Time-dependent solution  

 

 

> T[n](t):=C(n)*exp(k*lambda[n]*t);
 

`*`(C(n), `*`(exp(`+`(`-`(`*`(`/`(1, 40), `*`(`^`(n, 2), `*`(`^`(Pi, 2), `*`(t)))))))))
 

>
 

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

`*`(C(n), `*`(exp(`+`(`-`(`*`(`/`(1, 40), `*`(`^`(n, 2), `*`(`^`(Pi, 2), `*`(t))))))), `*`(cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x)))))))))
 

> u[0](x,t):=C(0)*X[0](x);
 

`+`(`*`(`/`(1, 2), `*`(C(0), `*`(`^`(2, `/`(1, 2))))))
 

>
 

Eigenfunction expansion  

 

 

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

`+`(`*`(`/`(1, 2), `*`(C(0), `*`(`^`(2, `/`(1, 2))))), Sum(`*`(C(n), `*`(exp(`+`(`-`(`*`(`/`(1, 40), `*`(`^`(n, 2), `*`(`^`(Pi, 2), `*`(t))))))), `*`(cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x)))))))...
 

 

Evaluation of Fourier coefficients for the specific initial condition  

 

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

`*`(x, `*`(`+`(2, `-`(x))))
 

 

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

 

Int(`*`(x, `*`(`+`(2, `-`(x)), `*`(cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x))))))))), x = 0 .. 2)
 

 

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

 

`+`(`-`(`/`(`*`(8, `*`(`+`(1, `^`(-1, n)))), `*`(`^`(n, 2), `*`(`^`(Pi, 2))))))
 

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

Int(`+`(`*`(`/`(1, 2), `*`(x, `*`(`+`(2, `-`(x)), `*`(`^`(2, `/`(1, 2))))))), x = 0 .. 2)
 

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

`+`(`*`(`/`(2, 3), `*`(`^`(2, `/`(1, 2)))))
 

>
 

>
 

>
 

Check Fourier coefficients  

 

 

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

Plot_2d
 

>
 

 

Our simple solutions 

 

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

`+`(`-`(`/`(`*`(8, `*`(`+`(1, `^`(-1, n)), `*`(exp(`+`(`-`(`*`(`/`(1, 40), `*`(`^`(n, 2), `*`(`^`(Pi, 2), `*`(t))))))), `*`(cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x)))))))))), `*`(`^`(n, 2), `*`(`^...
 

> u[0](x,t):=C(0)*X[0](x);
 

`/`(2, 3)
 

 

Series solution 

 

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

 

`+`(`/`(2, 3), Sum(`+`(`-`(`/`(`*`(8, `*`(`+`(1, `^`(-1, n)), `*`(exp(`+`(`-`(`*`(`/`(1, 40), `*`(`^`(n, 2), `*`(`^`(Pi, 2), `*`(t))))))), `*`(cos(`+`(`*`(`/`(1, 2), `*`(n, `*`(Pi, `*`(x)))))))))), `*...
 

 

First few terms of sum 

 

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

 

`+`(`/`(2, 3), `-`(`/`(`*`(4, `*`(exp(`+`(`-`(`*`(`/`(1, 10), `*`(`^`(Pi, 2), `*`(t)))))), `*`(cos(`*`(Pi, `*`(x)))))), `*`(`^`(Pi, 2)))))
 

Animation 

 

> animate(plot,[u(x,t),x=0..L],t=0..5,color=red,thickness=3, numpoints=200, frames=30);
 

 

 

 

Plot_2d
 

The preceding animation command illustrates the spatial-time dependent solution for u(x, t). The animation sequence  shows snapshots at times t = 0, 1, 2, 3, 4, and 5. 

 

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);
 

Plot_2d
 

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

Plot_2d
 

>
 

 

Example 3 Heat Equation Dirichlet-Neumann Conditions 

 

Assignment of system parameters  

 

> restart:
 

> L:=1; k:=1/4;
 

 

1
`/`(1, 4)
 

> with(plots):
 

>
 

Eigenvalues and orthonormal eigenfunctions.  

 

 

 

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

 

`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(n)), `-`(1)), `*`(Pi))))
`+`(`-`(`*`(`/`(1, 4), `*`(`^`(`+`(`*`(2, `*`(n)), `-`(1)), 2), `*`(`^`(Pi, 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)):
 

 

`*`(`^`(2, `/`(1, 2)), `*`(sin(`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(n)), `-`(1)), `*`(Pi, `*`(x))))))))
 

>
 

Statement of orthonormality  

 

 

 

 

 

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

 

Int(`+`(`*`(2, `*`(sin(`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(n)), `-`(1)), `*`(Pi, `*`(x)))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(m)), `-`(1)), `*`(Pi, `*`(x)))))))))), x = 0 .. 1) = delta...
 

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

1
 

Time-dependent solution  

 

 

> T[n](t):=C(n)*exp(k*lambda[n]*t);
 

`*`(C(n), `*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(`+`(`*`(2, `*`(n)), `-`(1)), 2), `*`(`^`(Pi, 2), `*`(t)))))))))
 

>
 

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

`*`(C(n), `*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(`+`(`*`(2, `*`(n)), `-`(1)), 2), `*`(`^`(Pi, 2), `*`(t))))))), `*`(`^`(2, `/`(1, 2)), `*`(sin(`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(n)), `-`(1)), `*`(...
 

>
 

>
 

Eigenfunction expansion  

 

 

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

Sum(`*`(C(n), `*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(`+`(`*`(2, `*`(n)), `-`(1)), 2), `*`(`^`(Pi, 2), `*`(t))))))), `*`(`^`(2, `/`(1, 2)), `*`(sin(`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(n)), `-`(1)), ...
 

 

Evaluation of Fourier coefficients for the specific initial condition  

 

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

`*`(x, `*`(sin(`*`(Pi, `*`(x)))))
 

 

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

 

Int(`*`(x, `*`(sin(`*`(Pi, `*`(x))), `*`(`^`(2, `/`(1, 2)), `*`(sin(`+`(`*`(`/`(1, 2), `*`(`+`(`*`(2, `*`(n)), `-`(1)), `*`(Pi, `*`(x)))))))))), x = 0 .. 1)
 

 

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

 

`+`(`/`(`*`(4, `*`(`^`(2, `/`(1, 2)), `*`(`+`(`-`(`*`(8, `*`(n))), 4, `-`(`*`(3, `*`(`^`(-1, n), `*`(Pi)))), `*`(4, `*`(`^`(-1, n), `*`(`^`(n, 2), `*`(Pi)))), `-`(`*`(4, `*`(`^`(-1, n), `*`(n, `*`(Pi)...
 

>
 

Check Fourier coefficients  

 

 

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

Plot_2d
 

>
 

 

Our simple solutions 

 

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

`+`(`/`(`*`(8, `*`(`+`(`-`(`*`(8, `*`(n))), 4, `-`(`*`(3, `*`(`^`(-1, n), `*`(Pi)))), `*`(4, `*`(`^`(-1, n), `*`(`^`(n, 2), `*`(Pi)))), `-`(`*`(4, `*`(`^`(-1, n), `*`(n, `*`(Pi)))))), `*`(exp(`+`(`-`(...
`+`(`/`(`*`(8, `*`(`+`(`-`(`*`(8, `*`(n))), 4, `-`(`*`(3, `*`(`^`(-1, n), `*`(Pi)))), `*`(4, `*`(`^`(-1, n), `*`(`^`(n, 2), `*`(Pi)))), `-`(`*`(4, `*`(`^`(-1, n), `*`(n, `*`(Pi)))))), `*`(exp(`+`(`-`(...
 

 

Series solution 

 

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

 

Sum(`+`(`/`(`*`(8, `*`(`+`(`-`(`*`(8, `*`(n))), 4, `-`(`*`(3, `*`(`^`(-1, n), `*`(Pi)))), `*`(4, `*`(`^`(-1, n), `*`(`^`(n, 2), `*`(Pi)))), `-`(`*`(4, `*`(`^`(-1, n), `*`(n, `*`(Pi)))))), `*`(exp(`+`(...
Sum(`+`(`/`(`*`(8, `*`(`+`(`-`(`*`(8, `*`(n))), 4, `-`(`*`(3, `*`(`^`(-1, n), `*`(Pi)))), `*`(4, `*`(`^`(-1, n), `*`(`^`(n, 2), `*`(Pi)))), `-`(`*`(4, `*`(`^`(-1, n), `*`(n, `*`(Pi)))))), `*`(exp(`+`(...
 

 

First few terms of sum 

 

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

 

`+`(`/`(`*`(`/`(8, 9), `*`(`+`(`-`(4), `*`(3, `*`(Pi))), `*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(Pi, 2), `*`(t)))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(Pi, `*`(x))))))))), `*`(`^`(Pi, 2))), `/`(`*`(`/`(8...
`+`(`/`(`*`(`/`(8, 9), `*`(`+`(`-`(4), `*`(3, `*`(Pi))), `*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(Pi, 2), `*`(t)))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(Pi, `*`(x))))))))), `*`(`^`(Pi, 2))), `/`(`*`(`/`(8...
`+`(`/`(`*`(`/`(8, 9), `*`(`+`(`-`(4), `*`(3, `*`(Pi))), `*`(exp(`+`(`-`(`*`(`/`(1, 16), `*`(`^`(Pi, 2), `*`(t)))))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(Pi, `*`(x))))))))), `*`(`^`(Pi, 2))), `/`(`*`(`/`(8...
 

Animation 

 

> animate(plot,[u(x,t),x=0..L],t=0..10,color=red,thickness=3, numpoints=200, frames=30);
 

 

 

 

Plot_2d
 

>
 

The preceding animation command illustrates the spatial-time dependent solution for u(x, t). The animation sequence  shows snapshots at times t = 0, 1, 2, 3, 4, and 5. 

 

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);
 

Plot_2d
 

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

Plot_2d
 

>