Mathematica Examples for Math 330

Heat Equation

Define the heat equation:

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

Initial condition: a box function

Define the initial condition:

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

Solve the differential equation using DSolve:

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

Plot the solution:

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

Initial condition: a parabola; homogeneous Dirichlet boundary conditions

Define the initial conditions:

ic=u[x,0]== x-x^2;

Define the boundary conditions:

bc1=u[0,t]==0;
bc2=u[1,t]==0;

Solve the differential equation using NDSolve:

sol=NDSolve[{heatEq,ic,bc1,bc2},u,{x,0,1},{t,0,1}]

Plot the solution:

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

Wave Equation

Run the following to clear all variables:

Clear["Global`*"]

Define the wave equation:

waveEq=D[u[x,t],{t,2}]==D[u[x,t],{x,2}]

Initial condition: a smooth bump; no boundary conditions (infinite string)

Define the initial conditions:

ic={u[x,0]==Exp[-x^2],Derivative[0,1][u][x,0]==0}

Solve the differential equation using DSolve:

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

Plot the solution:

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

Now impose boundary conditions: we will fix endpoints at zero (homogeneous Dirichlet boundary conditions)

Define the boundary conditions:

bc={u[-5,t]==0,u[5,t]==0}

Solve the differential equation using NDSolve:

sol=NDSolve[{waveEq,u[x,0]==Exp[-x^2],Derivative[0,1][u][x,0]==0,u[-5,t]==0,u[5,t]== 0},u,{x,-5,5},{t,0,20}]

Plot the solution. The following plot is colored by height.

Plot3D[Evaluate[u[x,t]/.sol[[1]] ],{x,-5,5},{t,0,20},PlotRange->All, AxesLabel->All,ColorFunction->Hue]

Wave equation on a circular disk

Define the circularly-symmetric wave equation:

circWaveEq=r D[u[r,t],{t,2}]==D[r D[u[r,t],r],r];

Initial conditions:

ic={u[r,0]==1,Derivative[0,1][u][r,0]==r/3};

Boundary conditions:

bc=u[1,t]==0;

Solve the wave equation using DSolve:

sol=DSolve[{circWaveEq,ic,bc},u[r,t],{r,t}]

Note that the solution is given as an infinite sum of Bessel functions.

Rewrite the solution:

sol2=(u[r,t]/.sol // FullSimplify)

Select the first three terms of the sum:

h[r_,t_]=sol2/.{Infinity-> 5} //Activate//N

Plot the first three terms of the sum at fixed time t:

Plot3D[ h[r,0.7]/.{r-> Sqrt[x^2+y^2]},{x,y}∈Disk[],PlotRange->All]

Animate[Plot3D[ h[r,t]/.{r-> Sqrt[x^2+y^2]},{x,y}∈Disk[],PlotRange->{-5,5}],{t,0,5}]

Poisson equation on a rectangular region

From the Wolfram tutorial Solving Partial Differential Equations with Finite Elements.

Needs["NDSolve`FEM`"]

Example 1

sol=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==10,u[x,0]==0,u[x,1]==1},u,{x,0,1},{y,0,1}]

Plot the solution:

Plot3D[sol[x,y],{x,0,1},{y,0,1},ColorFunction->"TemperatureMap",AxesLabel->Automatic]

We can also display the mesh on which the numerical solution was computed:

Show[ContourPlot[sol[x,y],{x,0,1},{y,0,1},AspectRatio->Automatic],sol["ElementMesh"]["Wireframe"]]

Example 2

For this example, we will define a region on which we will solve the PDE:

Ω=ImplicitRegion[True,{{x,0,2},{y,0,1}}]

Now define boundary conditions:

Γ=DirichletCondition[u[x,y]==0,x==0||x==2]

Solve the PDE:

sol=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==20,Γ},u,{x,y}∈Ω]

Plot the solution:

Plot3D[sol[x,y],{x,y}∈Ω,ColorFunction->"TemperatureMap",AxesLabel->Automatic]

Here is a contour plot that shows the mesh used by NDSolveValue:

Show[ContourPlot[sol[x,y],{x,y}∈Ω,AspectRatio->Automatic],sol["ElementMesh"]["Wireframe"]]

Poisson equation on an irregular region

From the Wolfram tutorial Solving Partial Differential Equations with Finite Elements.

First, define the region:

Ω=ImplicitRegion[!((x-5)^2+(y-5)^2≤3^2),{{x,0,5},{y,0,10}}];
RegionPlot[Ω,AspectRatio->Automatic]

Boundary conditions:

Subscript[Γ,D]={DirichletCondition[u[x,y]==0,x==0&&8≤y≤10],DirichletCondition[u[x,y]==100,(x-5)^2+(y-5)^2==3^2]};

Solve the PDE:

sol=NDSolveValue[{-Laplacian[u[x,y],{x,y}]==20,Subscript[Γ,D]},u,{x,y}∈Ω]

Visualize the solution:

ContourPlot[sol[x,y],{x,y}∈Ω,ColorFunction->"Temperature",AspectRatio->Automatic,PlotLegends->Automatic]

Plot3D[ufun[x,y],{x,y}∈Ω,ColorFunction->"Temperature",AspectRatio->Automatic]