MATH529

SP.24: Homework 4

Instructions

This homework assignment considers extension of the idea of separation of variables to more complicated, non-homogeneous partial differential equations.

Previously established results can be used directly in the solution.

Exercises

There are no exercises in this assignment

Problems

From p.746:6-16

Model solution:  13.6.5 (p. 746)

Statement

Solve the IBVP: 
(Subscript[u, t])=k (Subscript[u, x x])+ A (e^(-β x)), for 0<x<1, t>0, A∈ℝ
u(x,t=0)=f(x) for 0<=x<=1
u(x=0,t)=0 for t>0
u(x=1,t)=0 for

Solution

This is a non-homogeneous heat IBVP with a source term A (e^(-β x)), and homogeneous Dirichlet boundary conditions. Attempt a solution by superposition u(x,t)=v(x,t)+w(x) with associated problems

vIBVP: (Subscript[v, t])=k (Subscript[v, x x]), v(0,t)=0, v(1,t)=0, v(x,0)=g(x)=f(x)-w(x)

wBVP: k w''+A (e^(-β x))=0, w(0)=0, w(1)=0

Integrating the ODE in wBVP twice leads to

w'=(A/k β) (e^(-β x)) +(Subscript[c, 1]), w=-(A/ k (β^2)) (e^(-β x))+(Subscript[c, 1]) x+(Subscript[c, 2])

Applying boundary  conditions:

w(0)=-A/ k (β^2)+ (Subscript[c, 2])=0⇒(Subscript[c, 2])=A/ k (β^2),  
w(1)=-(A/ k (β^2)) (e^(-β))+(Subscript[c, 1])+A/ k (β^2)=0⇒(Subscript[c, 1])=A/ k (β^2)((e^(-β))-1),  hence

w(x)=(A/ k (β^2)) [1-(e^(-β x))-x(1-(e^(-β)))]

The vIBVP is homogeneous with homogeneous Dirichlet boundary conditions, with a sine series solution

v(x,t)=(∑_(n=1)^∞)(Subscript[B, n]) sin(n π x) (e^(-(((n π))^2) t/k)), and coefficients determined by

(Subscript[B, n])=2 (∫_0^1)g(x)sin(n π x)dx, g(x)=f(x)-(A/ k (β^2)) [1-(e^(-β x))-x(1-(e^(-β)))].

The overall solution is

u(x,t)=(∑_(n=1)^∞)(Subscript[B, n]) sin(n π x)(e^(-(((n π))^2) t/k))+(A/ k (β^2)) [1-(e^(-β x))-x(1-(e^(-β)))]

Solution study

As mentioned, solution verification and visualization is informative, not required, but does aid in mastery of course concepts.

Solution verification:

Clear["Global`*"];
wODE = k w''[x]+ A (E^(-β x))==0; wBC0 = w[0]==0; wBC1 = w[1] ==0; wBVP = {wODE,wBC0,wBC1}

wsol=Simplify[DSolve[wBVP,w[x],x][[1,1]]]

uPDE = (∂_t)u[x,t]== k (∂_(x,x))u[x,t]+A (E^(-β x))

u[x_,t_]= B Sin[n π x] (E^(-(((n π))^2) k t))+w[x] /. wsol

Assuming[n∈Integers,Simplify[uPDE]]

v[x_,t_,nt_]:= Sum[B[n] Sin[n π x] (E^(-(((n π))^2) k t)),{n,1,nt}];

u[x_,t_,nt_]:=v[x,t,nt]+w[x] /. wsol;

u[x,t,3]

f[x_]=x(1-x); g[x_]=f[x] - w[x] /. wsol

Bcoef[n_]:=Assuming[ n∈Integers, Integrate[g[x] Sin[n π x],{x,0,1}] ]

Bcoefs=Table[Bcoef[n],{n,1,5}]

B[n_]:=Bcoefs[[n]]

u[x_,t_,A_,k_,β_]=u[x,t,5]

Heating the bar

Plot3D[u[x,t,1,0.1,0.1],{x,0,1},{t,0,10},AxesLabel->{"x","t","u"}]

Cooling the bar

Plot3D[u[x,t,-0.1,0.1,0.1],{x,0,1},{t,0,10},AxesLabel->{"x","t","u"}]

Projects

From p.750: 6-10. Include a visualization and comments of your solution

Model solution: 13.7.2 (p. 750)

Statement

Solve the IBVP: 
(Subscript[u, t])=k (Subscript[u, x x]), for 0<x<1, t>0,
u(x,t=0)=f(x) for 0<=x<=1
u(x=0,t)=0 for t>0
(Subscript[u, x])(1,t)=-h(u(1,t)-(Subscript[u, 0]))⇔(Subscript[u, x])(1,t)+h u(1,t) = h (Subscript[u, 0]) for

Solution

This is a heat IBVP with fixed in time boundary conditions, homogeneous Dirichlet at left and nonhomogeneous Robin at right. The Robin BC states that the heat flux at x=1 is proportional to the difference of the bar and ambient medium temperatures.

Construct a solution by superposition u(x,t)=v(x,t)+w(x) (time-independent Robin condition). Note: this is the essential step in the solution, recognition of time-independent boundary conditions.

For w(x) formulate wBVP :
(Subscript[w, x x])=0, for 0<x<1, w(0)=0,  (Subscript[w, x])(1)+h w(1) =h (Subscript[u, 0])

with solution

w(x)= h (Subscript[u, 0]) x/(1+h)

For v(x,t) formulate an IBVP with homogeneous BC’s, vIBVP:

(Subscript[v, t])=k (Subscript[v, x x]), for 0<x<1, t>0,
v(x,0)=g(x)=f(x)-w(x) for 0<=x<=1
v(0,t)=0 for t>0
(Subscript[v, x])(1,t)+h v(1,t) =0 for 

Attempt a separation of variables solution u(x,t)=X(x)T(t) and obtain

(1/k)((T')/T)=((X'')/X)=λ

For λ=0, X(x)=(Subscript[c, 1])+(Subscript[c, 2])x, X(0)=0⇒(Subscript[c, 1])=0,
X'(1)+ h X(1)=(Subscript[c, 2])+h (Subscript[c, 2])=0⇒(Subscript[c, 2])=0, null solution.

For λ=(α^2)>0, solution becomes unbounded as t→∞, non-physical.

For λ=-(α^2)<0, X(x)=(Subscript[c, 3]) cos(α x)+(Subscript[c, 4]) sin(α x). BC’s imply:
X(0)=0=(Subscript[c, 3]), X'(1)+h X(1)=0=(Subscript[c, 4]) [α cos α+h sin α]⇒
α cos α+h sin α=0⇒tan α=-α/h, an eigencondition, with solutions (Subscript[α, n]), n∈ℕ\{0}. Construct solution as 

v(x,t)=(∑_(n=1)^∞)(Subscript[A, n]) (e^(- ((Subscript[α, n])^2) k t)) sin ((Subscript[α, n]) x)

with coefficients (Subscript[A, n]) determined to satisfy the initial condition

v(x,0)=(∑_(n=1)^∞)(Subscript[A, n]) sin((Subscript[α, n]) x)=g(x).

Overall solution is

u(x,t)=(∑_(n=1)^∞)(Subscript[A, n]) (e^(- ((Subscript[α, n])^2) k t)) sin ((Subscript[α, n]) x)+h (Subscript[u, 0]) x/(1+h).

Solution visualization

Arbitrarily choose some (Subscript[A, n]) coefficients, e.g., 
(Subscript[A, n])=1/(n^2) (this is equivalent to arbitrarily choosing some f(x), but simpler)

Clear["Global`*"];
A[n_] = 1/(n^2);

The eigenvalues can be found by numerically solving tan α + α/h=0. Choose h=1, and plot tan α+α

Plot[Tan[α]+α,{α,0,20}]

A first non-zero root is around α≃2

FindRoot[Tan[α]+α==0,{α,2}]

The next root is around 5

FindRoot[Tan[α]+α==0,{α,5}]

Subsequent roots are expected to be found at intervals of π, the period of tan. Construct a set of eigenvalues.

αEV = Table[ α /. FindRoot[Tan[α]+α==0,{α,4.9+ k π}],{k,-1,20}]

nT = Length[αEV]

Construct the solution

u[x_,t_,k_,u0_] = Sum[ A[n] Exp[-(((αEV[[n]]))^2)  k t] Sin[ αEV[[n]] x] + u0 x/2,{n,1,nT}];

Plot3D[u[x,t,1,0.1],{x,0,1},{t,0,1},AxesLabel->{"x","t","u"}]

Comment: The initial condition has higher temperatures than (Subscript[u, 0]). After a transient the temperature stabilizes to a linear function between endpoints. This is seen in time cross-sections of above plot.

Plot[Table[u[x,t,1,0.1],{t,0.,1.,0.2}],{x,0,1}]

Additional discussion

The eigencondition tan α=α/h does not lead to a Fourier series (known to be orthogonal). Orthogonality can be checked numerically to avoid tedious trigonometric identities

scProd[m_,n_]:=NIntegrate[Sin[αEV[[m]] x] Sin[αEV[[n]] x],{x,0,1}]

P=Table[scProd[m,n],{m,1,nT},{n,1,nT}];

The matrix of scalar products has non-negligible terms only on the diagonal.

ArrayPlot[P]

Note however that the basis is orthogonal, but not orthonormal

nrms=Diagonal[P]

The above scaling factors have to be taken into account when finding the (Subscript[A, n]) coefficients for some specific initial condition f(x).

f[x_]=x(1-x)

Af = Table[ NIntegrate[f[x] Sin[αEV[[n]]x ]/nrms[[n]],{x,0,1}],{n,1,nT}]

fSeries[x_] =Sum[ Af[[n]] Sin[αEV[[n]] x],{n,1,nT}];

Plot[{f[x],fSeries[x]},{x,0,1}]