Parabolic Equations

George E. Hrabovsky

MAST

Introduction

In this lesson we look, in some detail, at the properties of parabolic equations. In specific we will study the heat and diffusion equations. We will explore many of the same issues we covered in the last lesson and then we will examine convection-type problems. We will examine systems with losses. We will also cover solutions over specific regions. We will end with a discussion of the time-dependent Schrödinger equation.

Heat Equation

Let’s say that we have a rod of length l. The temperature at any point along the rod is a function of the position along that rod, x, and time, t. Each material has its own ability to allow heat to flow through it, and this is represented by the symbol D (for diffusion, as we will see later). The equation that governs this flow of heat is called the heat equation, and is an iconic parabolic equation.

l10_1.png

(10.1)

We can establish Dirichlet conditions at the end points, by assigning l10_2.png as the fixed temperature at x=0 and l10_3.png for the fixed temperature at x=l.

l10_4.png

(10.2)

For any real material the rod will start at some initial temperature, l10_5.png, giving us the initial condition

l10_6.png

(10.3)

The nature of D varies depending upon the problem. For simple problems this is a parameter determined by the nature of the material. It might be a function of one or more variables. In some cases it is a tensor. While the heat equation in (10.1) might look fairly simple, the nature of D can make it arbitrarily complicated.

The Method of Separation of Variables

We can illustrate a method of solution called the separation of variables. If we combine the heat equation (10.1) with the boundary conditions of (10.2) and the initial conditions of (10.3) then we have what is called an Initial-Value Boundary-Value Problem (IVBVP). If we fix the temperature at the end points by the condition

l10_7.png

(10.4)

The solution of this problem is then of the form

l10_8.png

(10.5)

This solution is separable. If we substitute the solution into (10.1), we get,

l10_9.png

(10.6)

Dividing through by D f(x) g(t), we get

l10_10.png

(10.7)

These ratios are constant, we can assign the symbol λ to these, so

l10_11.png

(10.8)

This in turn leads us to two ODEs,

l10_12.png

(10.9)

and

l10_13.png

(10.10)

The solutions for these are given below.

l10_14.png

l10_15.png

and

l10_16.png

l10_17.png

If we consider (10.4) we can write this f(0)=0=f(l).

l10_18.png

l10_19.png

So l10_20.png. We can then find the constant for x=l.

l10_21.png

l10_22.png

l10_23.png

l10_24.png

Naively, we can solve this to get the constant.

l10_25.png

l10_26.png

The problem with this is that we are seeking a nontrivial solution. So we need to look at exp2. This implies that

l10_27.png

(10.11)

where n is some natural number greater than 0. We can substitute this into exp2 to get a general form for specific values of x.

l10_28.png

l10_29.png

Each of these values provides a solution to sol2.

We then have an infinite number of solutions based on (10.5), these solutions will be of the form

l10_30.png

(10.12)

By the principle of superposition we can write

l10_31.png

(10.13)

We can use this along with (10.3) to write

l10_32.png

(10.14)

So,

l10_33.png

(10.15)

are the eigenvalues for this problem. Similarly,

l10_34.png

(10.16)

are the eigenfunctions.

Say that we have two eigenfunctions whose integer values are not equal, then we have the result.

l10_35.png

l10_36.png

So these eigenfunctions are orthogonal. If the two eigenfunctions share the same integer order, we get this result.

l10_37.png

l10_38.png

From this we can write an expression for the constants,

l10_39.png

(10.17)

For any given eigenfunction we have this result.

l10_40.png

l10_41.png

We can insert this into (10.17).

l10_42.png

(10.18)

We can then substitute this into (10.13),

l10_43.png

(10.19)

We can check to see if we get the same answer using DSolve.

l10_44.png

l10_45.png

l10_46.png

l10_47.png

l10_48.png

l10_49.png

l10_50.png

This takes a long time. The result is clearly the same as (10.19).

Dirichlet Problem

So what happens if the Dirichlet conditions are not 0? We write,

l10_51.png

The solution for this is then given below.

l10_52.png

l10_53.png

Now we need to specify some values. Say that the rod has a length of 10 units, the initial temperature of the rod is l10_54.png with the endpoint temperatures are l10_55.png and l10_56.png. Let’s also give a heat diffusion parameter of D=1/20. So we write it this way.

l10_57.png

l10_58.png

We can then extract the first five terms of the series.

l10_59.png

l10_60.png

We can plot this.

l10_61.png

l10_62.gif

Steady State Problems

Over significant periods systems will achieve thermal equilibrium and the variation of temperature will fall away. So we can write,

l10_63.png

(10.20)

for some function of position. We can also write,

l10_64.png

(10.21)

We call this a steady state situation. The solution for this is then constructed.

l10_65.png

l10_66.png

We can rewrite this.

l10_67.png

l10_68.png

We can compare this to the limit of sol5 as t->,

l10_69.png

l10_70.png

l10_71.png

l10_72.png

Is this the same as exp10?

l10_73.png

l10_74.png

We can examine a plot of the previous situation for very long times.

l10_75.png

l10_76.gif

Diffusion

The heat equation is a special application of the diffusion equation,

l10_77.png

(10.22)

for some function of time and spatial coordinates, ρ, that represents density and where l10_78.png is a positive definite matrix of diffusion coefficients that may or may not depend on position and time—this is called the diffusion tensor.

We write the diffusion equation.

l10_79.png

l10_80.png

We establish some initial conditions.

l10_81.png

So can we solve this IVP?

l10_82.png

l10_83.png

This is a strange result, the density does not change through space and time. It is a bizarre prediction. Let’s try a numerical solution, given some specific value for l10_84.png. Also the diffusion coefficient will be treated as 1.

l10_85.png

l10_86.png

l10_87.png

l10_88.gif

The plot of this is given.

l10_89.png

l10_90.gif

Let’s try another initial condition along with a continuous source term.

l10_91.png

l10_92.png

l10_93.png

We can extract the first few terms of this series.

l10_94.png

l10_95.png

l10_96.png

l10_97.png

We can plot this.

l10_98.png

l10_99.gif

We can compare this to the numerical solution.

l10_100.png

l10_101.gif

l10_102.png

l10_103.gif

So, what happens if this medium is in the form of a ring, we impose periodic boundary conditions.

l10_104.png

We then get this solution.

l10_105.png

l10_106.png

We can extract the first few terms of this series.

l10_107.png

l10_108.png

l10_109.png

l10_110.png

We can plot this.

l10_111.png

l10_112.gif

We can compare this to the numerical solution.

l10_113.png

l10_114.png

l10_115.gif

l10_116.png

l10_117.gif

They are similar, but not the same. What happens if we take more terms in the series?

l10_118.png

l10_119.png

l10_120.gif

It has not changed, and it seems to almost violate the initial conditions along the x boundary, this is sensitive to the periodic boundary conditions. We can compare their difference.

l10_121.png

l10_122.gif

The problem is due to the lack of compliance with the boundary conditions in the numerical solution, these issues fade away as the model continues to evolve in time. We can attempt to use the default method for NDSolve, the method of lines, but we will impose a pseudospectral method for the difference order due to the periodicity of the solution.

l10_123.png

l10_124.png

l10_125.gif

The plot of this is then give here.

l10_126.png

l10_127.gif

Let’s see these two solutions plotted again st each other for t=0.

l10_128.png

l10_129.gif

This points out the discrepancy between a series type solution and a numerical solution. Are we seeing an inconsistency in the conditions of a numerical solution? Or are we seeing artifacts due to Gibbs phenomena? No way to tell, here.

Neumann Problems

Returning to the heat equation, we can impose the Neumann conditions.

l10_130.png

This gives us a new solution.

l10_131.png

l10_132.png

We can extract the first few terms of this series.

l10_133.png

l10_134.png

l10_135.png

l10_136.png

We can plot this.

l10_137.png

l10_138.gif

Inhomogeneous Equations

In general we write the inhomogeneous heat equation,

l10_139.png

(10.23)

We write this equation.

l10_140.png

If we use the initial and boundary conditions defined earlier we can construct the solution.

l10_141.png

l10_142.png

Note that the solution is in terms of additional integrals. One way to cope with these is to rewrite them as numerical integrals. We can specify some values and resolve the equation. Let’s assume that f(x,t)=sin x.

l10_143.png

l10_144.png

We can extract the first few terms of this series.

l10_145.png

l10_146.png

l10_147.png

l10_148.png

We can plot this.

l10_149.png

l10_150.gif

Convection

What happens if we have the same rod, but at the right end it experiences convective heat transfer. This convection is represented as a sink. We see this in the Robin-type boundary condition.

l10_151.gif

We can now write our solution.

l10_152.png

l10_153.png

This is another example of the method of separation of variables. The infinite series is over the eigenfunctions that correspond to the eigenvalues listed as K[2,K[1]]. These eigenvalues are defined by the condition,

l10_154.png

(10.24)

where K[2,K[1]]>0, K[1]≥1, and K[1]∈Z. If we want to plot our solution we need to calculate the eigenvalues that are within some range of values, say 0<K[2,K[1]]<500000.

l10_155.png

How many eigenvalues are there in this range?

l10_156.png

l10_157.png

So we need to extract the first 1125 terms of the series corresponding to the eigenvalues.

l10_158.png

l10_159.png

We need to get rid of the conditional statements at the end.

l10_160.png

l10_161.png

We now expand this by using the apply command @@.

l10_162.png

We can then plot this.

l10_163.png

l10_164.gif

Losses

How do we cope with the idea of a heat loss into the surrounding medium while fixing the temperature at both ends. If we represent this heat loss with the term η u, then we write,

l10_165.png

(10.25)

So we write the equation.

l10_166.png

We can write a set of boundary and initial values.

l10_167.gif

The solution is then given.

l10_168.png

l10_169.png

We can apply our values to this.

l10_170.png

l10_171.png

We can extract the first five terms of the series.

l10_172.png

l10_173.png

We can then plot this.

l10_174.png

l10_175.gif

Unbounded Domains

What happens if our mythical rod is extremely long? We can treat it as semi-infinite. While the length may tend to infinity, we assume that temperature remains finite. We will continue to use eqn1 for this model. We will then write the initial and boundary values.

l10_176.gif

We can then write the solution.

l10_177.png

l10_178.png

Let’s fill in some values,

l10_179.png

l10_180.png x>0
Indeterminate True

We can then plot this.

l10_181.png

l10_182.gif

l10_183.png

l10_184.gif

Two-Dimensions

What happens if we want to extend our study to a rectangular plate instead of a rod? It is the same kind of equation,

l10_185.png

(10.26)

To keep it simple let’s assume a square.

l10_186.png

We can establish our initial conditions.

l10_187.png

We have Dirichlet conditions for y=0 and y=10, with the temperature being zero there.

l10_188.png

We assume the x=0 has a Neumann boundary that is insulated with a temperature of 0, and the other end has a Robin condition at x=10 and is in contact with a low temperature surrounding medium.

l10_189.png

Instead of solving this symbolically and entering another rabbit hole of eigenvalues and eigenfunctions, I will solve this numerically.

l10_190.png

l10_191.png

l10_192.gif

We can plot this for any time.

l10_193.png

l10_194.gif

We can make an animation.

l10_195.png

l10_196.gif

Specific Regions

Here we define a region with a slit in it.

l10_197.png

What does this look like?

l10_198.png

l10_199.gif

We apply the same equation, eqn5. We will keep the initial condition.

l10_200.png

l10_201.png

We then construct a numerical solution.

l10_202.png

l10_203.png

l10_204.png

l10_205.png

So, we can rewrite our equation.

l10_206.png

l10_207.png

l10_208.gif

We test this for some specific time.

l10_209.png

l10_210.gif

We can make an animation of this.

l10_211.png

l10_212.gif

Time-Dependent Schrödinger Equation

We can write the Time-Dependent Schrödinger Equation (TDSE) in Mathematica.

l10_213.png

We will apply Dirichlet conditions for two endpoints, a,b.

l10_214.png

So, we can write the solution.

l10_215.png

l10_216.png

We can extract the first five terms.

l10_217.png

l10_218.png

As an exercise, determine if this solution satisfies the equation given any set of constants.

So, let’s specify the initial conditions, a=0,b=1 and we will apply the natural units for h, and mass will also be 1.

l10_219.png

l10_220.png

Created with the Wolfram Language