Learning PDEs with the Wolfram Language

Unal Goktas, Wolfram Research

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

An overview of solution methods for partial differential equations (PDEs) in the Wolfram Language and how these tools can be used to understand and learn PDEs will be presented. We will show examples of the methods in play as they are being used in the Wolfram Language, and give insight into the key ideas and algorithms supported by graphics illustrations. The talk will focus on the monograph for computing symbolic solutions of PDEs in Version 12.3 and on the new CompleteIntegral function in the upcoming release.

Outline

Simplest PDE

First Order Linear Equations-Method of Characteristics

Fixed End Temperatures and Separation of Variables

Steady-State Temperatures

Shock Waves and Rarefaction Fan Solutions

The Black-Scholes Equation

Spherical Harmonics

Modeling

New System Function: CompleteIntegral

Monge Cone

Conclusion

Simplest PDE

The simplest possible PDE is ((∂u)/(∂x))=0 (where u=u(x,y)). Its general solution is u=(Subscript[\[ConstantC], 1])(y), where (Subscript[\[ConstantC], 1]) is any function of one variable. For instance, u=y, u=(y^2), u=(y/(y+1)) and u=(e^y) are several solutions. Since the solutions do not depend on x, they are constant on the lines y=C in the xy plane, where C is a constant.

For this PDE, DSolve returns:

DSolve[(∂_({x}))u[x,y]==0,u[x,y],{x,y}]

The following Manipulate shows that the solutions are constant on the y=C curves:

Manipulate[Column[{Plot3D[(Subscript[f, 1]),{x,1,5},{y,1,1.05 C},],}],{(Subscript[f, 1]),{y,(y^2),(y/(y+1)),(E^y)}},{{C,2},2,6},]

First Order Linear Equations-Method of Characteristics

The general first-order linear PDE IVP with two independent variables is given as:

a(x,t) (Subscript[u, x])+b(x,t) (Subscript[u, t])+c(x,t)u =0, -∞<x<∞,

u(x,0)=ϕ(x), 0<t<∞.

One solution technique to solve first-order linear PDEs is the method of characteristics, where we aim to find a change of independent variables (x,t) to new variables (s,τ) in order to obtain an ODE IVP that is easier to solve than the original PDE IVP.

At this point we have to answer the question: How do we find the characteristic curves that reduce the PDE IVP into an ODE IVP?

In order to find the characteristic curve and solve the resulting ODE uniquely, we need to define an IVP. We require that new coordinates s (to parameterize x and t) and τ (for the initial curve) have the following properties:

s varies along the characteristic curves {(x(s),t(s)): 0<s<∞}, ranging from 0 to ∞ starting at the initial curve, the line t=0 for this case.

τ varies along the initial curve t=0 and stays constant along each characteristic curve.

The following figure illustrates these properties:

Hence to uniquely define the characteristic curves in the new variables s and τ, we define the IVPs:

((dx)/(ds))=a(x,t), x(0)=τ

and

((dt)/(ds))=b(x,t), t(0)=0.

Since t(0)=0, we have t(s)=0 when s=0. Hence in the new coordinate system, (x,0) is mapped onto (x(s=0),t(s=0))=(x(0),t(0))=(τ,0). Therefore, the initial condition u(x,0)=ϕ(x) is transformed into u(0)=ϕ(τ) to act as an initial condition for the reduced ODE in s.

By applying the chain rule and the preceding properties, we arrive at the ODE IVP:

((du)/(ds))+c(x(s),t(s)) u=0,

u(0)=ϕ(τ), 0<s<∞.

So along the characteristics curves {(x(s),t(s)): 0<s<∞}, the PDE problem is reduced down to an ODE IVP. Once we solve this newly constructed problem, we change back to the original independent variables (x,t) to obtain the solution to the original PDE.

In summary, this approach, where:

New variables (s,τ) define characteristic curves (x(s),t(s),u(s)), each of which has (τ,0) as the initial point, and the initial data propagates along these curves,

The s-coordinate indicates points along a characteristic curve,

The τ-coordinate indicates the initial point (t=0) for a given curve,

For each τ, the initial value ϕ(τ) is evolved along the curve starting at the point (τ,0) according to the ODE obtained by expressing the PDE in characteristic variables,

is known as the method of characteristics.

The characteristic equations, (Subscript[x, s])=a(x,t), (Subscript[t, s])=b(x,t), are autonomous. The term autonomous is used when the equation has no explicit dependence on s.

The characteristics satisfy the ODE:

((Subscript[t, s])/(Subscript[x, s]))=(((dt)/(ds))/((dx)/(ds)))=((dt)/(dx))=((b(x,t))/(a(x,t))).

Hence the characteristics can be solved independently of the solution u.

Note that if c≡0 in the PDE IVP, then u remains constant along the characteristic curves.

Example: First Order PDE with variable coefficients

The equation

(Subscript[u, x])+y (Subscript[u, y])=0

is linear and homogeneous but has a variable coefficient (y), for which the characteristic curves satisfy the ODE ((dy)/(dx))=(y/1)=y. Solving the ODE gives:

DSolve[y'[x]==y[x],y[x],x]

y=(e^x) C are the characteristics. u(x,y) is a constant on each such curve. So u(x,y)=f(C), where f is an arbitrary function. Therefore, the general solution is obtained by solving y=(e^x)C for C. That is:

Solve[y==(E^x) C,C]

Hence

u(x,y)=f((e^(-x)) y).

We get the same solution directly from DSolve:

DSolve[(∂_({x}))u[x,y]+y (∂_({y}))u[x,y]==0,u[x,y],{x,y},GeneratedParameters->((Subscript[f, #1])&)]

Geometrically, the "picture" of the solution u(x,y) is that it is constant on each characteristic curve:

Manipulate[Column[{Plot3D[(Subscript[f, 1]),{x,0,2},{y,1,12},],}],{(Subscript[f, 1]),{(E^(-x)) y,((((E^(-x)) y))^2),((((E^(-x)) y))^3)}},{{C,2},2,10},]

Find the solution that satisfies the auxiliary condition u(0,y)=(y^3). Indeed, putting x=0, we get (y^3)=f((e^(-0)) y), so that f(y)=(y^3). Therefore, u(x,y)=((((e^(-x)) y))^3)=(e^(-3 x)) (y^3). The same solution can be found by DSolve with the call:

DSolve[(∂_({x}))u[x,y]+y (∂_({y}))u[x,y]==0&&u[0,y]==(y^3),u[x,y],{x,y}]

Example: First-order PDE IVP with constant coefficients

Suppose we want to solve the IVP with constant coefficients:

(Subscript[u, x])+(Subscript[u, t])+2 u=0,

u(x,0)=sin(x), -∞<x<∞, 0<t<∞.

First, to find the characteristic curves, solve the following ODE IVPs:

((dx)/(ds))=1, x(0)=τ,

and

((dt)/(ds))=1, t(0)=0,

for which DSolve gives:

DSolve[(∂_s)x[s]==1&&x[0]==τ,x[s],s]

DSolve[(∂_s)t[s]==1&&t[0]==0,t[s],s]

Hence, our characteristic curves are:

x=s+τ

and

t=s.

Using this new coordinate system, we reduce the PDE to the ODE:

((du)/(ds))+2 u=0, u(0)=sin(τ).

Solving this IVP, we get:

DSolve[(∂_s)u[s]+2 u[s]==0&&u[0]==Sin[τ],u[s],s]

So

u(s)=(e^(-2 s)) sin(τ),

is the solution in the sτ-coordinate system. In order to find u as a function of x and t, solve the transformation

x=s+τ,

t=s,

for s and τ in terms of x and t. We have:

Solve[{x==s+τ,t==s},{s,τ}]

Hence, our answer is

u(x,t)=(e^(-2 t)) sin(x-t),

which we could have directly found in an equivalent form by DSolveValue:

DSolveValue[(∂_({x}))u[x,t]+(∂_({t}))u[x,t]+2 u[x,t]==0&&u[x,0]==Sin[x],u[x,t],{x,t}]

Plot3D[%,{x,-10,10},{t,0,3},]

If we plot the solution at three different time values, we get:

Plot[Evaluate[%%/.{{t->0},{t->1/2},{t->1}}],{x,-10,10},]

The preceding plot shows the initial wave u(x,0)=sin(x) in orange. Note that with time progressing to t=1/2 and t=1, the wave dampens out; in other words, the amplitude decreases. Also, note that the peaks at the various times move slightly to the right.

Fixed End Temperatures and Separation of Variables

Suppose that we have a rod with length l whose lateral sides are wrapped with insulation, so that heat can flow in and out of the rod at the ends, but not across the lateral boundary. If this rod is placed in an environment whose temperature is fixed at some temperature (Subscript[T, 0]), after a sufficient time, the temperature of the entire rod stabilizes and comes to a steady-state temperature similar to the environment. Then suppose we take the rod out of the environment at a time that we call t=0 and attach two temperature components to the ends of the rod to keep the ends at specific fixed temperatures, say (Subscript[T, 1]) and (Subscript[T, 2]). If we monitor the temperature chart of the rod, the basic PDE that is modeling the one-dimensional heat flow is:

(Subscript[u, t])=k (Subscript[u, xx]), 0<x<l, 0<t<∞.

PDE relates the quantities (Subscript[u, t]) and (Subscript[u, xx]), where the former represents the rate of change in temperature with respect to time, and the latter represents the concavity of the temperature profile u(x,t). In other words, PDE compares the temperature at one point to the temperature at neighboring points. The proportionality constant k is a property of the material that the rod is made of.

Since the temperature u was fixed for all times t>0 at (Subscript[T, 1]) and (Subscript[T, 2]) at the two endpoints x=0 and x=l, for the boundary conditions, we could say:

u(0,t)=(Subscript[T, 1])

and

u(l,t)=(Subscript[T, 2]),

i.e. the ends of the rod are held at fixed temperatures .

Recall that we started monitoring the rod temperature from the time the rod had reached a temperature of (Subscript[T, 0]); we have the initial condition:

u(x,0)=(Subscript[T, 0]), 0≤x≤l.

Hence the PDE, Dirichlet BCs and the IC make up the IVBVP.

Let us first suppose that (Subscript[T, 1])=(Subscript[T, 2])=0, i.e. the ends of the rod are held at zero temperature. So, we want to solve the heat flow IVBVP:

(Subscript[u, t])=k (Subscript[u, xx]), 0≤x≤l, 0<t<∞,

IC: u(x,0)=(Subscript[T, 0]),

BCs: u(0,t)=0 and u(l,t)=0.

We assume that the solution u(x,t) can be written as:

u(x,t)=X(x) T(t),

i.e. the solution is separable.

Substituting this form into the PDE leads to:

X(x) T'(t)=k T(t) X''(x).

Dividing both sides by k X(x) T(t) results in:

((T'(t))/(k T(t)))=((X''(x))/(X(x)))=-λ,

where λ has to be constant since X(x) is a function of solely x and T(t) is a function of solely t.

Now we have two ODEs:

T'(x)+k λ T(t)=0

and

X''(x)+λ X(x)=0,

for which the general solutions are:

T(t)=(e^(-k λ t))

and

X(x)=A cos((Sqrt[λ]) x)+B sin((Sqrt[λ]) x),

for λ>0, where A and B are constants.

The next step is imposing BCs, which requires that X(0)=0=X(l). Thus we have:

X(0)=A=0

and

X(l)=B sin((Sqrt[λ]) l)=0.

Since we are interested in a nontrivial solution, i.e. A and B should not both be zero, we arrive at:

sin((Sqrt[λ]) l)=0 ⇒(Sqrt[λ]) l=π n ⇒λ=(((((π n)/l)))^2), n=1,2,3,…

So for (Subscript[λ, n])=(((((π n)/l)))^2),

(Subscript[X, n])(x)=B sin((Sqrt[(Subscript[λ, n])]) x)=B sin(((n π x)/l)), n=1,2,3,…,

are distinct solutions.

Hence there are infinitely many solutions of the PDE such that:

(Subscript[u, n])(x,t)=(Subscript[B, n]) (e^(-k (((((n π)/l)))^2) t)) sin(((n π x)/l)), n=1,2,3,…,

where (Subscript[B, n]) are constants. Due to the superposition principle for linear equations, the sum of solutions is also a solution. Therefore,

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

is also a solution.

Imposing IC, we get:

(Subscript[T, 0])=(∑_(n=1)^∞)(Subscript[B, n]) sin(((n π x)/l)).

In the preceding, (Subscript[λ, n])=(((((π n)/l)))^2) are called eigenvalues and sin((Sqrt[(Subscript[λ, n])]) x)=sin(((n π x)/l)) are called eigenfunctions.

Using the orthogonality of the eigenfunctions, i.e.:

Assuming[(n|m)∈(Integers)&&n≥1&&m≥1&&n≠m,(∫_0^ℓ)Sin[((n π x)/ℓ)] Sin[((m π x)/ℓ)]\[DifferentialD]x]

Assuming[(n|m)∈(Integers)&&n≥1&&m≥1&&n==m,(∫_0^ℓ)Sin[((n π x)/ℓ)] Sin[((m π x)/ℓ)]\[DifferentialD]x]

and:

Assuming[n∈(Integers)&&n≥1,(∫_0^ℓ)Sin[((n π x)/ℓ)] \[DifferentialD]x]

we find that:

(Subscript[B, n])=(2/l)(Subscript[T, 0]) (∫_0^l)sin(((n π x)/l))dx=((2 (1+(((-1))^(n+1))) (Subscript[T, 0]))/(π n))=-((2 (-1+(((-1))^n)) (Subscript[T, 0]))/(π n)).

Substituting this into the form of the solution, we conclude that:

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

which is the same solution that DSolveValue returns:

pde=(∂_({t}))u[x,t]==k (∂_({x,2}))u[x,t];

ic=u[x,0]==(Subscript[T, 0]);

bcs=u[0,t]==0&&u[ℓ,t]==0;

solution=DSolveValue[pde&&ic&&bcs,u[x,t],{x,t}]

The method outlined is the so-called separation of variables method, also called the eigenfunction expansion method.

Note that the λ=0 and λ<0 cases do not lead to nontrivial solutions, and hence we omitted these cases.

Now let us return to the fixed-end temperatures problem. For this IVBVP, DSolveValue returns:

pde=(∂_({t}))u[x,t]==k (∂_({x,2}))u[x,t];

ic=u[x,0]==(Subscript[T, 0]);

bcs=u[0,t]==(Subscript[T, 1])&&u[ℓ,t]==(Subscript[T, 2]);

solution=DSolveValue[pde&&ic&&bcs,u[x,t],{x,t}]

For k=(1/10), ℓ=10, (Subscript[T, 0])=5, (Subscript[T, 1])=30 and (Subscript[T, 2])=60, a three-dimensional plot of the solution is:

k=(1/10);ℓ=10;(Subscript[T, 0])=5;(Subscript[T, 1])=30;(Subscript[T, 2])=60;

Plot3D[Evaluate[DSolveTruncatedSolution[solution,500]],{x,0,ℓ},{t,0,30},]

where, the function DSolveTruncatedSolution is defined as:

DSolveTruncatedSolution[solution_,truncateAt_Integer]:=Activate[solution/.∞->truncateAt]

Notice that while plotting the solution we truncated the infinite summation at K[1]=500 to get a sufficiently good approximation for the summation by using the function DSolveTruncatedSolution. Clearly, if we truncate at larger values of K[1], we could get better approximations. In this example, truncating at K[1]=500 was sufficient for our purpose.

In the next subsection, we will see that the first two terms in the preceding solution form the steady-state solution.

Steady-State Temperatures

Suppose that the function u(x,t) represents temperature. After a sufficiently long time under the same environment and conditions, experiences show that the variation of temperature with time dies away. Hence, we expect that the limit of u(x,t), as t tends to infinity, exists and depends only on x:

(Limit2Arg[u(x,t), t, ∞])=v(x)

and also that

(Limit2Arg[((∂u(x,t))/(∂t)), t, ∞])=0.

The steady-state temperature distribution function v(x) has to satisfy the BCs and the heat equation, which are valid for all t>0. Therefore, v(x) should be the solution to the problem:

DSolveValue[{v''[x]==0,v[0]==(Subscript[T, 1]),v[ℓ]==(Subscript[T, 2])},v[x],x]//(MapAll[Apart[#,x]&,#]&)

which is equal to the limit of the solution of the preceding IVBVP as t->∞.

Let k=(1/10), l=30, (Subscript[T, 0])=10, (Subscript[T, 1])=20 and (Subscript[T, 2])=40 in the preceding IVBVP and solve the problem with concrete values:

pde=(∂_({t}))u[x,t]==(1/10) (∂_({x,2}))u[x,t];

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

bcs=u[0,t]==20&&u[30,t]==40;

solutionfixedendtemperatures=DSolveValue[pde&&ic&&bcs,u[x,t],{x,t}]

The following Manipulate computation shows that after a long enough time, the solution starts converging to the steady-state solution:

Manipulate[Quiet[Show[
Plot3D[Evaluate[N[DSolveTruncatedSolution[solutionfixedendtemperatures,500]]],{x,0,30},{t,(Subscript[t, 0]),(Subscript[t, 0])+10},],
Graphics3D[{Opacity[.4],Blue,InfinitePlane[{{0,(Subscript[t, 0])+5,0},{1,(Subscript[t, 0])+5,0},{1,(Subscript[t, 0])+5,1}}]}]]],{{(Subscript[t, 0]),1000},1000,3000,50},]

Note that Quiet in the computation is used to suppress General::munfl messages that get generated during the computation.

Shock Waves and Rarefaction Fan Solutions

The IVP in which the initial data is formed of two constant states (Subscript[u, l]) and (Subscript[u, r]) and has a jump discontinuity at x=0 is called the Riemann problem. Now let us consider Burgers' equation with the Riemann problem:

(Subscript[u, t])+u (Subscript[u, x])=(Subscript[u, t])+(Subscript[((1/2)(u^2)), x])=0,

where the IC is given as:

u(x,0)=({{{, ({{(Subscript[u, l]), x<0}, {(Subscript[u, r]), x>0}})}}).

For (Subscript[u, l])>(Subscript[u, r]), the characteristics intersect, hence a smooth solution cannot be constructed. Instead, we can construct a weak solution of the form of a shock wave moving at a speed determined by the Rankine–Huguenot condition. For (Subscript[u, l])<(Subscript[u, r]), we can construct a weak solution of the form of a rarefaction fan (transonic rarefaction) from the following equation:

pde=(∂_t)u[x,t]+u[x,t] (∂_x)u[x,t]==0;

ic=u[x,0]==Piecewise[{{ul,x<0},{ur,x>0}}];

DSolveValue[pde&&ic,u[x,t],{x,t}]

So when u(x,0)=(UnitStepSeq[-x]), we have a shock wave solution:

pde=(∂_t)u[x,t]+u[x,t] (∂_x)u[x,t]==0;

ic=u[x,0]==UnitStep[-x];

solution=DSolveValue[pde&&ic,u[x,t],{x,t}]

for which a three-dimensional plot of the solution is:

Plot3D[solution,{x,-3,3},{t,0,5},]

Whereas when u(x,0)=(UnitStepSeq[x]), we have a rarefaction fan solution:

pde=(∂_t)u[x,t]+u[x,t] (∂_x)u[x,t]==0;

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

solution=DSolveValue[pde&&ic,u[x,t],{x,t}]

for which a three-dimensional plot of the solution is:

Plot3D[solution,{x,-3,3},{t,0,5},]

Example: Generalized Burgers’ Equation

Consider the generalized viscous Burgers' equation from the theory of turbulence:

(Subscript[u, t])+u (Subscript[u, x])=ϵ (Subscript[u, xx]),

where the IC is given as:

u(x,0)=({{{, ({{1, x<0}, {0, x⩾0}})}}).

For this IVBVP, DSolveValue returns:

pde=(∂_t)u[x,t]+u[x,t] (∂_x)u[x,t]==ϵ (∂_({x,2}))u[x,t];

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

solution=DSolveValue[pde&&ic,u[x,t],{x,t}]

The solution is smooth for any positive value of ϵ:

Plot3D[solution/.ϵ->(1/10),{x,-2,2},{t,0,5},]

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

Quiet[GraphicsRow[Table[Plot3D[solution,{x,-2,2},{t,0,5},],{ϵ,{(1/10),(1/100),(1/1000)}}]]]

The Black–Scholes Equation

All European options satisfy the Black–Scholes PDE when the underlying asset price is assumed to be risk-neutral. If V(s,t) is the value of the option with underlying asset s and elapsed time t, risk-free interest rate r, dividend yield q and volatility or annualized standard deviation of the asset price σ, then the equation modeling all European options is:

(Subscript[V, t])(s,t)+(r-q) s (Subscript[V, s])(s,t)+(1/2) (σ^2) (s^2) (Subscript[V, ss])(s,t)-r V(s,t)=0.

If k is the strike and T is the time of maturity, the IC for the put option is given as:

V(s,T)=max(k-s,0).

For the put option problem, DSolveValue returns:

pde=(∂_t)v[s,t]+(1/2) (σ^2) (s^2) (∂_({s,2}))v[s,t]+(r-q) s (∂_s)v[s,t]-r v[s,t]==0;

ic=v[s,T]==Max[k-s,0];

DSolveValue[pde&&ic,v[s,t],{s,t},Assumptions->σ>0&&0<q<r&&0<t<T&&k>0]

For the put option, the BCs are implicitly taken as:

V(0,t)=k (e^(-r (T-t)))

and

(Limit2Arg[V(s,t), s, ∞])=0.

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

%/. {t->0,s->100,k->100,σ->0.7,T->1,r->0.05,q->0.01}

Compare with the value given by FinancialDerivative:

FinancialDerivative[{"European","Put"},{"StrikePrice"->100.,"Expiration"->1},{"InterestRate"->0.05,"Dividend"->0.01,"Volatility"->0.7,"CurrentPrice"->100}]

The IC for the call option is given as:

V(s,T)=max(s-k,0).

For the call option problem, DSolveValue returns:

pde=(∂_t)v[s,t]+(1/2) (σ^2) (s^2) (∂_({s,2}))v[s,t]+(r-q) s (∂_s)v[s,t]-r v[s,t]==0;

ic=v[s,T]==Max[s-k,0];

DSolveValue[pde&&ic,v[s,t],{s,t},Assumptions->σ>0&&0<q<r&&0<t<T&&k>0]

For the call option, the BCs are implicitly taken as:

V(0,t)=0

and

(Limit2Arg[((V(s,t))/s), s, ∞])=1.

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

%/. {t->0,s->100,k->100,σ->0.7,T->1,r->0.05,q->0.01}

Compare with the value given by FinancialDerivative:

FinancialDerivative[{"European","Call"},{"StrikePrice"->100.,"Expiration"->1},{"InterestRate"->0.05,"Dividend"->0.01,"Volatility"->0.7,"CurrentPrice"->100}]

Spherical Harmonics

Suppose that we want to find the potential inside a sphere with radius 3 when the potential, g(θ,ϕ)=sin(θ) sin(ϕ), is given on the boundary. For this we have to solve the Laplace PDE-BVP in a sphere with no azimuthal symmetry:

pde=Laplacian[u[r,θ,ϕ],{r,θ,ϕ},"Spherical"]==0;

bc=u[3,θ,ϕ]==Sin[θ] Sin[ϕ];

solutionSphericalHarmonics=DSolveValue[pde&&bc,u[r,θ,ϕ],{r,θ,ϕ}]

Let us define the truncated solution:

Quiet[truncatedsolutionSphericalHarmonics[r_,θ_,ϕ_]=Chop[DSolveTruncatedSolution[solutionSphericalHarmonics[[1,1,1]],1]/. Integrate->NIntegrate]];

truncatedsolutionSphericalHarmonics[r,θ,ϕ]

Notice that truncating the summation at larger values of K[1] also gives the same solution:

Quiet[Chop[DSolveTruncatedSolution[solutionSphericalHarmonics[[1,1,1]],2]/. Integrate->NIntegrate]]

Quiet[Chop[DSolveTruncatedSolution[solutionSphericalHarmonics[[1,1,1]],3]/. Integrate->NIntegrate]]

Quiet[Chop[DSolveTruncatedSolution[solutionSphericalHarmonics[[1,1,1]],4]/. Integrate->NIntegrate]]

Here is a visualization of the spherical three-dimensional plot of the solution:

Manipulate[SphericalPlot3D[Chop[N[Abs[truncatedsolutionSphericalHarmonics[r,θ,ϕ]]]],{ϕ,0,2 π},{θ,0,π},],{r,0.0001,3},]

Modeling

Heat transfer is a discipline of thermal engineering that is interested in modeling of transmission of energy. Let us solve a slightly modified version of the basic heat transfer example:

The model domain Ω of width W=0.02` [m] and height H=0.01` [m] is a ceramic strip that is embedded in a high-thermal-conductive material. The side boundaries of the strip are maintained at a constant temperature (Subscript[T, hot])=1173 [K]. The top surface of the strip is losing heat via thermal convection, and thermal radiation to the ambient environment at (Subscript[T, amb])=323 [K] is omitted. The bottom boundary, however, is assumed to be thermally insulated.

Our goal is to find the steady-state temperature distribution of the ceramic strip.

Set up a rectangular domain with a width of W=0.02` [m] and a height of H=0.01` [m].

Ω2D=Rectangle[{0,0},{0.02,0.01}];

The thermal conductivity k, heat transfer coefficient h, density ρ, heat capacity (Subscript[C, p]) and emissivity ε of the ceramic trip are given by:

({{k, =, 3, [W/(m·K)]}, {h, =, 50, [W/((m^2)·K)]}, {ρ, =, 3690, [kg/(m^3)]}, {(Subscript[C, p]), =, 880, [J/(kg·K)]}, {ε, =, 0.7,  }})

Define the model parameters:

vars={T[x,y],{x,y}};

pars=<|"ThermalConductivity"->3,"AmbientTemperature"->323|>;

Set up temperature surface boundary conditions at the left and right boundaries:

(Subscript[Γ, temp])=HeatTemperatureCondition[x==0||x==0.02,vars,pars,<|"SurfaceTemperature"->1173|>]

Set up the convective boundary condition on the top surface:

(Subscript[Γ, convective])=HeatTransferValue[y==0.01,vars,pars,<|"HeatTransferCoefficient"->50|>]

A default thermally insulated boundary condition is implicitly applied on the remaining bottom boundary:

(Subscript[Γ, insulation])=NeumannValue[0,y==0]

Solve the heat transfer PDE model with NDSolveValue and DSolveValue:

pde={HeatTransferPDEComponent[vars,pars]==(Subscript[Γ, convective])+(Subscript[Γ, insulation]),(Subscript[Γ, temp])};

TfunNDSolve = NDSolveValue[pde, T, {x, y} ∈ Ω2D]

TfunDSolve=DSolveValue[pde,T,{x,y}∈Ω2D]//Chop

Visualize the steady-state temperature distribution:

ContourPlot[TfunNDSolve[x,y],{x,y}∈Ω2D,]

Quiet[ContourPlot[Evaluate[N[DSolveTruncatedSolution[TfunDSolve[x,y],10000]]],{x,y}∈Ω2D,]]

Plot3D[TfunNDSolve[x,y],{x,y}∈Ω2D,]

Quiet[Plot3D[Evaluate[N[DSolveTruncatedSolution[TfunDSolve[x,y],10000]]],{x,y}∈Ω2D,]]

CompleteIntegral

A complete integral of a first-order partial differential equation (PDE) in n variables is a solution that depends on n arbitrary constants (Subscript[c, 1]),(Subscript[c, 2]),…,(Subscript[c, n]).

A solution to the PDE that satisfy a specific initial condition can be obtained by constructing the envelope of a smoothly varying subfamily of simple solutions that depend on n-1 parameters, as illustrated below.

Example: Clairaut Equation

Find a complete integral of the Clairaut equation:

ClairautEquation = u[x, y] == x*D[u[x, y], {x}] + y*D[u[x, y], {y}] + 
     (1/2)*(D[u[x, y], {x}]^2 + D[u[x, y], {y}]^2);

The complete integral is given by a 2-parameter family of planes:

csol=u[x,y]/.CompleteIntegral[ClairautEquation, u, {x, y}][[1]]

Select a one-parameter family of these planes:

psol = csol/.{C[2]->-1}

Find the envelope of this one-parameter family of planes:

Solve[(D[psol,C[1]])==0,C[1]]

envelope=psol /.%[[1]]//Together

Verify that the envelope is also a solution:

ClairautEquation/.{u->Function[{x,y},Evaluate[envelope]]}//Simplify

Visualize the one-parameter family of planes and the envelope solution:

Animate[Plot3D[{(psol/.{C[1]->c}),envelope }, {x, -20, 20}, 
   {y, -20, 20}, ], {c, -6, 6}, 
  ]

Monge Cone

Suppose that we have a first-order partial differential order for w=w(x,y), of the form:

F(x,y,u,p,q)=0,

where p=(Subscript[u, x]) and q=(Subscript[u, y]).

A complete integral, being a relation between x, y and w, is the equation of a surface. Since it contains two arbitrary parameters, it belongs to a double infinite system of surface or to a single infinite system of family of surfaces.

CompleteIntegral[D[w[x,y],x]^2+D[w[x,y],y]^2==1&&w[0,0]==0,w[x,y],{x,y}]

completeintegralsolution1[x_,y_]=%[[1,1,2]];
completeintegralsolution2[x_,y_]=%%[[2,1,2]];

With[{body1=Quiet[Chop[N[completeintegralsolution1[x,y]]]],body2=Quiet[Chop[N[completeintegralsolution2[x,y]]]]},Manipulate[Show[Plot3D[body1,{x,-5,5},{y,-5,5},],Plot3D[body2,{x,-5,5},{y,-5,5},],],{(C[2]),-1,1},]]

All the tangent planes to possible integral surfaces through a point P form a family of planes enveloping a conical surface which is called Monge Cone or elementary cone of at the point P with P as its vertex . In other words, the problem of solving the original PDE is to find surfaces which touch the Monge cone at each point along a generator. Let us illustrate this property:

solutionCompleteIntegral[{x_,y_},{c1_,c2_}]=c1+y c2-x (Sqrt[1-(c2^2)]);

tangentPlane[{x_,y_},{x0_,y0_},{c1_,c2_}] =Solve[(ResourceFunction["TangentPlane"][solutionCompleteIntegral[{x,y},{c1,c2}],{x,x0},{y,y0},z,"StandardFormEquation"]//Simplify),z][[1,1,2]]

Now fix a w=(Subscript[w, 0]). This restricts the two parameters (Subscript[c, 1]), (Subscript[c, 2]) by implying a dependence:

solutionC1=Solve[solutionCompleteIntegral[{x0,y0},{c1,c2}] == w0, c1]

Compute the implicit form of the cone as the envelope of all tangent planes:

implicitTangentPlane[{x_, y_},{x0_,y0_,w0_}, c2_] =z-tangentPlane[{x,y},{x0,y0},{solutionC1[[1,1,2]],c2}]

coneSurface[{x_, y_, z_}, {x0_,y0_,w0_}]=GroebnerBasis[{implicitTangentPlane[{x, y},{x0,y0,w0}, c2],D[implicitTangentPlane[{x, y},{x0,y0,w0}, c2], c2],0},
{},{c2}][[1]]//Factor

We see that is a genuine quadratic.

Solution surface through point ((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])):

solutionZ[{x_, y_},{x0_,y0_,w0_},c2_]=tangentPlane[{x,y},{x0,y0},{solutionC1[[1,1,2]],c2}]

Manipulate[
Block[{coneSize=3,planeSize=4,x0y0w0Point,sPlot,cone,tPlane},
x0y0w0Point = Graphics3D[{Red, Sphere[{x0,y0,w0}, 0.2]}];
sPlot=Plot3D[Evaluate[solutionZ[{x, y},{x0,y0,w0},c2]], {x, -5,5},{y,-5,5},
];
cone = ContourPlot3D[Evaluate[coneSurface[{x, y, z}, {x0,y0,w0}]==0], {x,x0-coneSize,x0+coneSize},{y,y0-coneSize,y0+coneSize},{z,w0-coneSize,w0+coneSize},
];
tPlane = ContourPlot3D[Evaluate[implicitTangentPlane[{x, y},{x0,y0,w0}, c2]==0], {x,x0-planeSize,x0+planeSize},{y,y0-planeSize,y0+planeSize},{z,w0-planeSize,w0+planeSize},
];
Show[{sPlot,x0y0w0Point,cone,tPlane},]
],
{{x0, 1}, -5,5},{{y0, 1}, -5,5},{{w0, 1}, -5,5},
{{c2, 0.3}, -1,1},]

Conclusion

For details and for more examples, see the manuscript: