Symbolic PDEs and Eigenfunctions
Devendra Kapadia
Wolfram Research, Inc.

Join the Conversation #WolframTechConf

Outline of the Talk

I would like to discuss the following new features of Version 10.3 during this talk:

Solutions of classical PDEs - Wave, Heat, and Laplace

Solutions of modern PDEs - Burgers’, Black-Scholes, Schrödinger

Sturm-Liouville Problems

Eigensystems of Differential Operators

Wave Equation

Initial Value Problem

The one-dimensional wave equation is:

WaveEquation=(((∂^2)u(x,t))/(∂(t^2)))==(((∂^2)u(x,t))/(∂(x^2)));

Prescribe initial conditions:

InitialCondition={u(x,0)==(E^(-(x^2))),Derivative[0,1][u](x,0)==1};

Solve the initial value problem:

sol=DSolve[{WaveEquation,InitialCondition},u[x,t],{x,t}]

Plot the solution:

Plot3D[Evaluate[u[x, t] /. %[[1]]], {x, -7, 7}, {t, 0, 4}, PlotRange -> All]

The wave equation is a basic example of a hyperbolic PDE.

A Vibrating String

Start with the wave equation:

eqn = D[u[x, t], {t, 2}] == D[u[x, t], {x, 2}];

Specify that the ends of the string remain fixed during the vibrations:

bc={u[0,t]==0,u[π,t]==0};

Give the initial position of the string:

ic={u[x,0]==(x^2) (π-x),Derivative[0,1][u][x,0]==0};

Solve the initial-boundary value problem:

dsol=DSolve[{eqn,bc,ic},u,{x,t}]/.{K[1]-> m}

Extract four terms from the Inactive sum:

asol[x_,t_]=u[x,t]/.dsol[[1]]/.{∞->4}//Activate

Each term in the sum represents a standing wave:

Table[Show[Plot[Table[asol[x,t][[m]],{t,0,4}]//Evaluate,{x,0,Pi},Ticks->False],ImageSize-> 150],{m,4}]

Visualize the vibration of the string:

Animate[Plot[asol[x,t],{x,0,π},PerformanceGoal->"Speed",PlotRange->{-5,5},ImageSize->Medium,PlotStyle->Red],{t,0,10},SaveDefinitions->True,AnimationRepetitions->3]

Heat Equation

Initial Value Problem

The one-dimensional heat equation is:

HeatEquation=((∂u(x,t))/(∂t))==(((∂^2)u(x,t))/(∂(x^2)));

Prescribe an initial condition:

InitialCondition=u[x, 0] == UnitBox[x];

Solve the initial value problem:

sol=DSolve[{HeatEquation,InitialCondition},u[x,t],{x,t}]

Plot the solution:

Plot3D[Evaluate[u[x,t]/.sol[[1]]],{x,-2,2},{t,0,1},PlotRange->All]

The heat equation is a basic example of a parabolic PDE.

Flow of Heat in a Bar

Model the flow of heat in a bar of length 1 that is insulated at both ends:

heqn=D[u[x,t],t]==D[u[x,t],{x,2}];

Specify that no heat flow through the ends of the bar:

bc={Derivative[1,0][u][0,t]==0,Derivative[1,0][u][1,t]==0};

Specify the initial temperature at different points in the bar:

ic=u[x,0]==20+80x;

Solve the heat equation subject to these conditions:

sol=DSolve[{heqn,bc,ic},u[x,t],{x,t}]

Extract a few terms from the Inactive sum:

approxsol=u[x,t]/.sol[[1]]/.{∞->6}//Activate//Expand

Visualize the evolution of the temperature to the steady-state value of 60 degrees:

Plot[Table[approxsol,{t,0,0.9,0.07}]//Evaluate,{x,0,1},AxesOrigin->{0,0},PlotRange->All]

Laplace Equation

The Dirichlet Problem in a Rectangle

The Laplace equation in two dimensions is:

LaplaceEquation=(Laplacian[u(x,y), {x,y}])==0;

Prescribe boundary conditions in a rectangle:

BoundaryConditions={u[x,0]==UnitTriangle[2 x-1],u[x,2]==UnitTriangle[2 x-1],u[0,y]==0,u[1,y]==0};

Solve the Dirichlet problem:

sol=(DSolve[{LaplaceEquation,BoundaryConditions},u[x,y], {x,y}][[1]]//FullSimplify)

Extract 300 terms from the Inactive sum:

asol=u[x,y]/.sol/.∞->300//Activate;

Visualize the solution:

Plot3D[asol,{x,0,1},{y,0,2},PlotRange->All,PlotStyle->Hue[0.8]]

The Laplace equation is a basic example of an elliptic PDE.

The Region Notation

Dirichlet problem for the Laplace equation in a rectangle using the region notation:

leqn=(Laplacian[u[x,y], {x,y}])==0;

dcond=DirichletCondition[u[x,y]==Piecewise[{{UnitTriangle[2 x-1],y==0||y==2}},0],True];

Ω=Rectangle[{0,0},{1,2}];

sol=FullSimplify[DSolve[{leqn,dcond},u[x,y],{x,y}∈Ω]]

Extract the first 300 terms from the Inactive sum:

asol=u[x,y]/.sol[[1]]/.{∞->300}//Activate;

Visualize the solution on the rectangle:

Plot3D[asol//Evaluate,{x,y}∈Ω,PlotRange->All,PlotStyle->Hue[0.8]]

Three Modern PDEs

Burgers’ Equation

Burgers' equation for viscous fluid flow is:

BurgersEquation=((∂u(x,t))/(∂t))+u(x,t) ((∂u(x,t))/(∂x))==ϵ (((∂^2)u(x,t))/(∂(x^2)));

Prescribe a piecewise initial condition:

InitialCondition=u[x,0]==Piecewise[{{1,x<0}}];

Solve the initial value problem:

dsol=DSolveValue[{BurgersEquation,InitialCondition}, u[x,t],{x,t}]

The solution is smooth for any positive value of ϵ:

Plot3D[dsol/.{ϵ-> 1/10},{x,-2,2},{t,0.001,5}]

The solution develops a shock discontinuity in the limit when ϵ approaches 0:

Table[Plot3D[dsol,{x,-2,2},{t,0.001,5},Exclusions->All,Ticks->None],{ϵ,{1/10,1/100,1/1000}}]

Black-Scholes Equation

The following is the classic Black-Scholes model from finance:

BlackScholesModel={((∂c(t,s))/(∂t))+(1/2) (s^2) (σ^2) (((∂^2)c(t,s))/(∂(s^2)))+r s ((∂c(t,s))/(∂s))-r c(t,s)==0,c(T,s)==max(s-k,0)};

Obtain the Black-Scholes formula:

dsol=DSolveValue[BlackScholesModel,c[t,s],{t,s}]

Compute the value of an option for typical values of the parameters:

dsol/.{t-> 0,s->100,k-> 100, σ-> 0.2,T-> 1,r-> 0.05}

Compare with the value given by FinancialDerivative:

FinancialDerivative[{"European","Call"}, {"StrikePrice"-> 100.00, "Expiration"->1},  {"InterestRate"-> 0.05, "Volatility" -> 0.2 ,"CurrentPrice"-> 100}]

Schrödinger Equation

Solve the free Schrödinger equation for an electron in a box:

SchrodingerEquation=I ℏ ((∂ψ(x,t))/(∂t))==-((ℏ^2)/(2 m)) (((∂^2)ψ(x,t))/(∂(x^2)));

BCS= {ψ[0,t]==0,ψ[d,t]==0};

ICS=ψ[x,0]==(1/(Sqrt[d]))Sin[Pi(x/d)]+(1/(Sqrt[d]))Sin[2Pi(x/d)];

sol2=DSolve[{SchrodingerEquation,BCS,ICS},ψ[x,t],{x,t}]

Compute the resulting probability density:

ρ2[x_,t_]=Simplify[ComplexExpand[Conjugate[ψ[x,t]]ψ[x,t] /. First[sol2]],d>0]

Substitute suitable parameter values:

ρe[x_,t_]=ρ2[x,t]/. {ℏ->1.055*^-34,m->9.109*10^-31,d->10^-9}

The "center" of the electron moves from side to side of the box:

Animate[Plot[ρe[x,t],{x,0,10^-9},PerformanceGoal->-"Quality",PlotRange->{0,4*^9},ImageSize->Medium,AxesLabel->{"x (m)", "ρe ((m^(-1)))"},Ticks->False],{t,0,3.66*^-15},SaveDefinitions->True,AnimationRepetitions-> 3]

Sturm-Liouville Problems

A Basic Eigenvalue Problem

Solve an eigenvalue problem with Dirichlet conditions:

sol=DSolve[{y''[x]+λ y[x]==0,y[0]==0,y[π]==0},y[x],x]

Make a table of the first 5 eigenfunctions:

eigfuns=Table[y[x]/.sol[[1]]//.{\[FormalN]->i,λ->(\[FormalN]^2)}/.{C[1]-> 1},{i,5}]

Plot the eigenfunctions:

Plot[Evaluate[eigfuns],{x,0,π}]

Transcendental Roots and Assumptions

Solve an eigenvalue problem with a Robin condition at the left end of the interval:

sol=DSolve[{y''[x]+λ^2 y[x]==0,y[0]+y'[0]==0,y[1]==0},y[x],x]

Locate the roots of the transcendental eigenvalue equation in the range 1<λ<10:

roots=λ/.Solve[λ Cos[λ]==Sin[λ]&&1<λ<10,λ];data=Table[{λ,0},{λ,roots}];

Plot[λ Cos[λ]-Sin[λ],{λ,1,10},Epilog->{PointSize[0.03],Red,Point[data]}]

Obtain the eigenfunctions within this range by using the new Assumptions option:

DSolveValue[{y''[x]+λ^2 y[x]==0,y[0]+y'[0]==0,y[1]==0},y[x],x,Assumptions->1<λ<10]

Eigensystems of Differential Operators

DEigensystem and DEigenvalues

DEigensystem and DEigenvalues can be used to compute eigenfunctions and eigenvalues of ordinary and partial differential operators.

For example:

DEigensystem[{-Laplacian[u[x],{x}],DirichletCondition[u[x]==0,True]},u[x],{x,0,π},4]

Visualize the eigenfunctions:

Plot[Evaluate[%[[2]]],{x,0,π}]

Obtain the eigenvalues only:

DEigenvalues[{-Laplacian[u[x],{x}],DirichletCondition[u[x]==0,True]},u[x],{x,0,π},4]

Eigenfunctions of the Laplacian in a Ball

Specify a Laplacian operator:

ℒ = -Laplacian[u[x, y, z], {x, y, z}];

Specify homogeneous Dirichlet boundary conditions :

ℬ=DirichletCondition[u[x,y,z]==0,True];

Find the 7 smallest eigenvalues and eigenfunctions in a ball:

{vals, funs} = DEigensystem[{ℒ, ℬ}, u[x, y, z], {x, y, z} ∈ Ball[{0, 0, 0}, 2], 7];

Visualize an eigenfunction:

DensityPlot3D[funs[[7]] // N // Evaluate, {x, y, z} ∈ Ball[{0, 0, 0}, 2], Boxed -> False, Axes -> False, ColorFunction -> Hue,ImageSize->Medium]

Conclusion

Version 10.3 introduces a new framework for symbolic PDEs and eigenfunctions.

These features rely on recent developments in visualization, numerics and symbolic computation.

The new functionality addresses a substantial number of requests made by our users.