MATH529

SP.22: Homework 6

Instructions

In this homework assignment you investigate how separation of variables applied to PDEs formulated in curvilinear coordinates leads to eigenfunctions different from the trigonometric and hyperbolic functions.

The ODEs that arise in such problem are difficult to solve. Historically, professional mathematicians have found solutions and introduced special functions for the widely encountered problems in orthogonal coordinate systems. Prior to the advent of symbolic software systems, such solutions would be found by table lookup. In this assignment, feel free to use DSolve to find solutions to Sturm-Liouville problems.

Exercises

There are no exercises in this assignment

Problems

From pp.760-1:11-20, p.768-9:14-18, p.771: 5-11

Model solution:  14.2.13 (p. 768)

Statement

Solve the IBVP: 
(Subscript[u, t])=(Subscript[u, r r])+(1/r)(Subscript[u, r]), for 0<r<2, t>0,
u(r,t=0)=f(r), 0<r<2 (initial condition) 
u(r=2,t)=100, t>0, (boundary condition)

f(r)=200 for 0<r<1, f(r)=100, for 1<r<2.

Solution

This is a heat IBVP  with non-homogeneous, constant in time Dirichlet boundary conditions.  The problem is formulated in cylindrical coordinates (with no variation along θ,z directions). Physical considerations exclude  infinite values for u(r,t). Attempt a solution by superposition u(r,t)=v(r,t)+w(r), with homogeneous Dirichlet conditions for v.

vIBVP: (Subscript[v, t])= (Subscript[v, r r])+(1/r) (Subscript[v, r]), v(r,0)=f(r)-w(r), v(2,t)=0

wBVP:  w''+(1/r) w'=0, w(2)=100

The wBVP contains a separable ODE. Denote p=w', and integrate p'+(1/r)p=0,
((dp)/(dr))=-(1/r) p⇒((dp)/p)=-((dr)/r)⇒ln p=-ln r+ln a⇒p=w'=a/r⇒w(r)=a ln r+b

At r=0, w finite implies a=0, while at r=2, w(2)=100, hence w(r)=100, a constant function.

Attempt a solution of the vIBVP by separation of variables

v(r,t)=R(r) T(t)⇒((T')/T)=((R'')/R)+(1/r) ((R')/R)

The separation constant must be negative to avoid infinite temperature values for t→∞, hence
T(t)∼(e^(-(α^2)t)), and the radial variation is described by ODE

r R''+R'+(α^2) r R=0,

a Sturm-Liouville problem, with a denumerable, complete set of eigenfunctions.

RODE = r R''[r]+R'[r]+(α^2) r R[r]==0; RBC0 = R[0]==0; RBC1 = R[2] ==0; RBVP = {RODE,RBC0,RBC1}

DSolve only returns the trivial null solution (is not coded to search for eigenvalues).

Rsol=Simplify[DSolve[RBVP,R[r],r][[1,1]]]

The general solution is in terms of Bessel functions

Rsol=Simplify[DSolve[RODE,R[r],r][[1,1]]]

To avoid infinite values at r=0, where (Subscript[J, 0])(0)=1, (Subscript[Y, 0])(0)=-∞, impose (Subscript[\[ConstantC], 2])=0, hence R(r)=(Subscript[\[ConstantC], 1]) (Subscript[J, 0])(α r). At r=2 the boundary condition v(2,t)=0 implies R(2)=(Subscript[\[ConstantC], 1]) (Subscript[J, 0])(2α)=0, which is satisfied non-trivially by choosing 2(Subscript[α, n]) to be roots of the Bessel function.

Plot[BesselJ[0,2r],{r,0,20} ,Frame->True,FrameLabel->{"2r","(Subscript[J, 0])(2r)"},GridLines->Automatic,PlotLabel->"Bessel function (Subscript[J, 0])(r)"]

α=Table[α /. FindRoot[BesselJ[0,2α],{α,3+1.5k}],{k,-1,20}]

The solution to the vIBVP is of form

v(r,t)=(∑_(k=1)^∞)(Subscript[A, k]) (Subscript[J, 0])((Subscript[α, k]) r) (e^(-((Subscript[α, k])^2) t))

The eigenfunction set is {(Subscript[J, 0])((Subscript[α, 1]) r),(Subscript[J, 0])((Subscript[α, 2]) r),...} and is orthogonal with respect to the scalar product

(f,g)=(∫_0^2)r f(r) g(r)dr

This can be verified for the first few eigenfunctions

S=Table[Integrate[r BesselJ[0,α[[j]] r] BesselJ[0,α[[k]] r],{r,0,2}],{j,1,5},{k,1,5}];
ArrayPlot[S]

The norms of the eigenfunctions can be computed

d=Monitor[Table[Integrate[r BesselJ[0,α[[k]] r] BesselJ[0,α[[k]] r],{r,0,2}],{k,1,22}],k]

The coefficients (Subscript[A, k]) are now evaluated from the initial condition

f[r_]:=If[r<1,200,100];
Plot[f[r],{r,0,2}]

A = Monitor[Table[Integrate[r BesselJ[0,α[[k]] r] (f[r]-100),{r,0,2}]/d[[k]],{k,1,22}],k];

The overall solution is

u(r,t)=v(r,t)+w(r)=(∑_(k=1)^∞)(Subscript[A, k]) (Subscript[J, 0])((Subscript[α, k]) r) (e^(-((Subscript[α, k])^2) t)) + 100

Solution study

Define the solution

u[r_,t_,nt_]:=Sum[A[[k]] BesselJ[0,α[[k]] r] Exp[-(((α[[k]]))^2) t],{k,1,nt}]+100

The initial condition discontinuity leads to Gibbs phenomenon and a solution that is not physically valid.

Plot[Table[u[r,t,5],{t,0.01,0.1,0.01}],{r,0,2} ,Frame->True,FrameLabel->{"r","u(r)"},GridLines->Automatic,PlotLabel->"Solution u(r,t) at t=0,0.01,...,0.1",PlotRange->All]

The temperature should never decrease below 100 or go over 200 in this problem. Increasing the number of terms does not help much

Plot[Table[u[r,t,22],{t,0,0.1,0.01}],{r,0,2} ,Frame->True,FrameLabel->{"r","u(r)"},GridLines->Automatic,PlotLabel->"Solution u(r,t) at t=0,0.01,...,0.1",PlotRange->All]

The series would converge to the initial condition as k→∞, not a particularly useful result. Try another, smooth initial condition to verify a correct procedure

f[r_]=100 (Cos[π r/4]+1);
Plot[f[r],{r,0,2}]

A = Monitor[Table[NIntegrate[r BesselJ[0,α[[k]] r] (f[r]-100),{r,0,2}]/d[[k]],{k,1,22}],k];

u[r_,t_,nt_]:=Sum[A[[k]] BesselJ[0,α[[k]] r] Exp[-(((α[[k]]))^2) t],{k,1,nt}]+100

With no initial discontinuity the solution is now physically correct and quite accurate with only a few terms in the series.

Plot[Table[u[r,t,5],{t,0,1,0.1}],{r,0,2} ,Frame->True,FrameLabel->{"r","u(r)"},GridLines->Automatic,PlotLabel->"Solution u(r,t) at t=0,0.1,...,0.1",PlotRange->All]

Projects

From pp.768-9: 20-22. Include a visualization and comments of your solution. No model solution is presented. Use approaches from Problems and HW05 Project to guide your study.