Solving PDEs with initial and boundary conditions:

Sturm-Liouville problem with RootOf eigenvalues

Computer algebra systems always failed to compute exact solutions for a linear PDE with initial / boundary conditions when the eigenvalues of the corresponding Sturm-Liouville problem cannot be solved exactly - that is, when they can only be represented at most abstractly, using a RootOf construction.

This post illustrates then a new Maple development, to tackle this kind of problem (work in collaboration with Katherina von Bülow), including testing and plotting the resulting exact solution. To reproduce the computation below in Maple 2018.1 you need to install the Maplesoft Physics Updates (version 134 or higher).

As an example, consider

This problem represents the temperature distribution in a thin rod whose lateral surface is insulated over the interval , with the left end of the rod insulated and the right end experiencing a convection heat loss, as explained in George A. Articolo's Partial Differential Equations and Boundary Value Problems with Maple, example 3.6.4.

The formulation as a Sturm-Liouville problem starts with solving the PDE by separating the variables by product

Substituting this separation by product into the last two (out of three) initial/boundary conditions (iv), the original pde and these conditions are transformed into an ODE system with initial conditions

This is a problem in actually three variables, including , and the solution can be computed using dsolve

We are interested in a solution of the form , so discard the first two in the above and keep only the third one

If we were able to express  in closed form, with a formula depending on an integer parameter identifying one of the roots of the expression, the rest of the method is straightforward, the product  is a solution that by construction satisfies the boundary conditions in (1) , and we have one of them for each value of (for each root of the expression within the RootOf construction). The first condition in iv is used to adjust the remaining constant (combine _C4 times _C5 into one) using orthogonality properties , and by linearly superimposing these different solutions we construct an infinite series solution. All this is explained in standard textbooks.

The problem is what to do when  cannot be expressed in closed form, as in the example above. To understand the situation clearly, remove that RootOf and plot its contents:

This equation has infinitely many roots. For instance between  and ,

A closed form solution to represent these values of  (also called the eigenvalue of the Sturm-Liouville problem) are necessary in the intermediate solving steps, and to express in closed form the solution to the original problem.

We cannot do magic to overcome this mathematical limitation, but there are in Maple representational abstract alternatives: we can use, in all the intermediate computations, a RootOf construction with a label  identifying each of the roots and, at the end, remove the RootOf construction, presenting something readable, understandable.

This is the representation of the solution that we are proposing, whose automatic derivation from given pde and iv is already implemented:

So an expression restricted by equations, inequations and possibly subject to conditions. This is the same type of representation we use in pdsolve, PDEtools:-casesplit and the FunctionAdvisor.

Note that, since there are infinitely many roots to the left and to the right of the origin, we assumed for simplicity that , which suffices to construct the infinite series solution above. The initial value for the summation index n could be 0 or 1, it doesn't matter, since n itself does not appear in the summand, it only identifies one of the values of  solving . This is a very nice development.

So we can now compute and represent the solution algebraically even when the eigenvalue  cannot be expressed in closed form. But how do you test such a solution? Or even plot it?

For now that needs some human intervention. Start with testing (part of what follows is in the plans for further development of pdetest). Separate the solution expressed in terms of from the equation defining it

By inspection, solution has  and , not , so rewrite (10), and on the way isolate

Start by pdetesting

A further manipulation, substituting condition and combining the sums results in

So from the four elements (one PDE and three iv), three are satisfied without having to specify who is  - it sufficed to say that , and this is the case of most of the examples we analyzed. You really don't need the exact closed form of  in these examples.

For the one expression which remains to be proven equal to zero, there is no clear way to perform the sum and show that it is equal to  without further information on the value of .  So this part needs to be tested using a plot.

We need to perform the summation, instead of using infinite terms, using, say, the first 100 of them, and for that purpose we need the first 100 positive values of . These values can be obtained using fsolve. Increase Digits to get a better approximation:

For convenience, construct a procedure, as a function of n, that returns each of these values

Replace  by infinity by 99 and the expression to be plotted is

Perform the sum and plot. The plotting range is the one present in the iv (1), x goes from 0 to 1

This result is very good. With the first 100 terms (constructed using the first 100 roots of ) we verified that this last of the four testing conditions is sufficiently close to zero, and this concludes the verification.

To plot the solution, the idea is the same one: in (9), give a value to k - say k = 1/5 - then construct the sum of the first 100 terms as done in (17), but this time using (9) instead of (13)

Compare with the numerical solution one could obtain using pdsolve with the numeric option . So substitute , and from the corresponding help page rewrite the initial/boundary conditions using D notation and this is the syntax and corresponding plot