• Example 1
  • Problem
  • Solving problem
  • Solution
• Example 2
  • Problem
  • Solving problem
  • Solution
• Example 3
  • Problem
  • Solving problem
  • Solution
• Example 4
  • Problem
  • Solving problem
  • Solution
  • Checking the Solution

Fourier Approximate Solutions to PDE Boundary Value Problems

by Aleksas Domarkas

Vilnius University, Faculty of Mathematics and Informatics,

Naugarduko 24, Vilnius, Lithuania

aleksas@ieva.mif.vu.lt

In this session we find  approximate solutions to boundary value problems for heat  and wave equations using Fourier method. We solve problems from    Numerical Solutions to PDE Boundary Value Problems in Maple 8.

Example 1

Problem

>	restart; Order:=30:

>	Heat_Eqn:=diff(u(x,t),t)-a^2*diff(u(x,t),x,x)=0; #0<x<1, t>0

>	BCs:={u(x,0)=phi(x),u(0,x)=0, u(1,x)=0};

>	a:=sqrt(1/10);

>	phi:=x->(1-x)*(1-exp(-50*x));

>	plot(phi,0..1);

>

Solving problem

We solve corresponding eigenvalue problem:

>	L:=-diff(y(x),x,x)=lambda*y(x);

>	s1:=y(0)=0; s2:=y(1)=0;

>	assume(lambda>0);

>	sol:=rhs(dsolve(L,y(x)));

>	sist:={subs(x=0,sol),subs(x=1,sol)};

>	linalg[genmatrix](sist,{_C1,_C2});

>	chl:=combine(linalg[det](%));

>	_EnvAllSolutions := true:

>	solve(chl,lambda);

>	indets(%);

>	subs(%[1]='k',%%):

>	ev:=unapply(%,k);

>	assume(k,posint);

>	subs(lambda=ev(k),L);

>	dsolve({%,s1,s2},y(x));

>	rhs(%)/sqrt(int(rhs(%)^2,x=0..1));

>	simplify(%,radical,symbolic):

>	ef:=unapply(%,(k,x)):

Eigenvalues and eigenfunctions:

>	ev(k);ef(k,x);

Solution we find in form:

>	solution:=Sum(T[k](t)*ef(k,x),k=1..Order);

For coefficients      we solve ODE problem:

>	Tk:=simplify(dsolve({diff(T[k](t),t)+a^2*ev(k)*T[k](t)=0,

>	T[k](0)=int(phi(x)*ef(k,x),x=0..1)}));

>	map(simplify,value(subs(Tk,solution))):

>

Solution

>	u:=unapply(%,(x,t)):

>	`u(x,t)`=sort(map(factor,u(x,t)),[seq(sin(Pi*k*x),k=1..Order)]);

The solution at  x=0.5 , t=3

>	evalf(u(0.5,3),20);

>	plot3d(u(x,t),x=0..1,t=0..2,axes=boxed,labels=["x","t","u(x,t)"]);

>	plots[animate](u(x,t),x=0..1,t=0..2,frames=30,labels=["x","u(x,t)"]);

>	plot([phi(x),u(x,0)],x=0..1);

>

Example 2

Problem

>

restart;

>	Heat_Eqn:=diff(u(x,t),t)-exp(-2*t)*rho(x)*diff(u(x,t),x,x)=0; #0<x<1, t>0

>	BCs:={u(x,0)=phi(x),u(0,x)=0, u(1,x)=0};

>	rho:=x->x/2+1/1000;

>	phi:=x->sin(Pi*x);

>

Solving problem

We solve corresponding eigenvalue problem:

>	L:=-diff(y(x),x,x)=lambda*y(x)/rho(x);

>	s1:=y(0)=0; s2:=y(1)=0;

>	assume(lambda>0);

>	sol:=rhs(dsolve(L,y(x)));

>	sist:={subs(x=0,sol),subs(x=1,sol)};

>	linalg[genmatrix](sist,{_C1,_C2});

>	combine(linalg[det](%));

>	ch_eq:=unapply(simplify(%/sqrt(12525)),lambda);

>	plot(ch_eq(x),x=1..50);

>	#plot(ch_eq(x),x=50..150);

The ranges  of roots  ch_eq(x)=0  are:

>	S:=1..5,5..10,10..15,20..30,30..40,40..60; #60..80,80..100,100..120,120..140;

>	m:=nops([S]);

We find eigenvalues and eigenfunctions:

>	for k to m do

>	mu[k]:=fsolve(ch_eq(x), x=S[k]); subs(_C2=1,lambda=mu[k],x=0,sol): isolate(%,_C1): subs(%,_C2=1,lambda=mu[k],sol): ef[k]:=evalf(%/sqrt(evalf(int(%^2/rho(x),x=0..1)))); print(lambda[k]=mu[k]);print(y[k]=ef[k]); end do: k:='k':

Solution we find in form:

>	solution:=Sum(T[k](t)*ef[k],k=1..m);

For coefficients      we solve ODE problems:

>	for k to m do

>	s[k]:=simplify(dsolve({diff(T[k](t),t)+exp(-2*t)*mu[k]*T[k](t)=0,

>	T[k](0)=evalf(int(phi(x)*ef[k]/rho(x),x=0..1))})); end do;

>

k:='k':

>

Solution

>	u:=unapply(subs(seq(s[k],k=1..m),value(solution)),(x,t));

The solution at x=0.5 , t=3

>	evalf(u(0.5,3),20);

>	plot3d(u(x,t),x=0..1,t=0..2,axes=boxed,labels=["x","t","u(x,t)"]);

>	plots[animate](u(x,t),x=0..1,t=0..2,frames=30,labels=["x","u(x,t)"]);

>

Example 3

Problem

>	restart; Order:=20:

>	Wave_Eqn:=diff(u(x,t),t,t)-a^2*diff(u(x,t),x,x)=0; #0<x<1, t>0

>	BCs:={u(x,0)=phi(x), D[2](u)(x,t)=0, u(0,x)=0, u(1,x)=0};

>	a:=sqrt(1/10);

case 1:

>	phi:=x->piecewise(x<=1/2,x,1-x);

>	plot(phi,0..1);

>

Solving problem

Corresponding egenvalue problem:

>	L:=-diff(y(x),x,x)=lambda*y(x);

>	s1:=y(0)=0; s2:=y(1)=0;

Eigenvalues and eigenfunctions(same as in Example 1):

>	ev:=unapply(Pi^2*k^2,k);

>	ef:=unapply(sin(Pi*k*x)*sqrt(2),(k,x));

Solution of this problem is sought in in the form:

>	solution:=Sum(T[k](t)*ef(k,x),k=1..Order);

For coefficients   T[k](t)  we solve ODE problem:

>	Tk:=simplify(dsolve({diff(T[k](t),t,t)+a^2*ev(k)*T[k](t)=0,

>	T[k](0)=int(phi(x)*ef(k,x),x=0..1),D(T[k])(0)=0}));

>	map(simplify,value(subs(Tk,solution))):

>

Solution

>	u:=unapply(%,(x,t));

>	plot3d(u(x,t),x=0..1,t=0..4*Pi,axes=boxed,labels=["x","t","u(x,t)"]);

>	plots[animate](u(x,t),x=0..1,t=0..2*Pi,frames=30,labels=["x","u(x,t)"]);

>	plot([phi(x),u(x,0)],x=0..1);

>

Example 4

Problem

>

restart;

>	Wave_Eqn:=diff(u(x,t),t,t)-a^2*diff(u(x,t),x,x)=0; #0<x<1, t>0

>	BCs:={u(x,0)=phi(x), D[2](u)(x,t)=0, u(0,t)=0, u(1,t)=0};

>	a:=sqrt(1/10);

case 2:

>	phi:=x->sin(4*Pi*x)/2;

>	plot(phi,0..1);

>

Solving problem

Corresponding eigenvalue problem:

>	L:=-diff(y(x),x,x)=lambda*y(x);

>	s1:=y(0)=0; s2:=y(1)=0;

Eigenvalues and eigenfunctions(same as in Example 1):

>	ev:=unapply(Pi^2*k^2,k);

>	ef:=unapply(sin(Pi*k*x)*sqrt(2),(k,x));

Solution of this problem is sought in in the form:

>	solution:=sum(T[k](t)*ef(k,x),k=1..Order);

For coefficients      we solve ODE problems:

>	for k to Order do

>	s[k]:=combine(dsolve({diff(T[k](t),t,t)+a^2*ev(k)*T[k](t)=0,

>	T[k](0)=int(phi(x)*ef(k,x),x=0..1),D(T[k])(0)=0}));od; k:='k':

>

Solution

In this case we have exact solution:

>	sol:=subs(seq(s[k],k=1..Order),solution);

>	u:=unapply(%,(x,t));

>	plot3d(u(x,t),x=0..1,t=0..Pi/2,axes=boxed,labels=["x","t","u(x,t)"]);

>	plots[animate](u(x,t),x=0..1,t=0..Pi/2,frames=30,labels=["x","u(x,t)"]);

>

Checking the Solution

>

Wave_Eqn;

>	u(x,0); is(%=phi(x));

>	D[2](u)(x,0);

>

u(0,t);

>

u(1,t);

>

>