Solving PDEs?



Saturday 23 June 2007 04:41:50 Thomas Widlar wrote:
> My friend has had good success with Maxima on his project. He makes the
> comment --
>
>   It seems Maxima is not able to handle partial differential equations
> (with several variables),
>   only ordinary differential equations. Mathematica and Maple have the
> reqired packages,
>   but are mainly tailored for mechanical and engineering applications.
>
>   Quantum mechanics is too special (requires algebraic eigenvalues for
> example). I think
>   if  feasable solutions of the Schr?dinger equation in cartesian
> coordinates existed, they
>   would be describe in books, at least in  the outline. I know that any
> group tried to solve
>   quantum mechanical equations numerically on a rectangular grid, but it
> was rather
>   unprecise because the grid was not adopted to the centralsymmetric
> geometry.
>
> Is he right about PDEs in Maxima?
That's right. 
Comercial macsyma has pdeasy. It's based on the finite elements methods and 
seems quite powerful. My personal opinion is that to create the general CAS 
interface to solve PDE numerically is a very hard task. Maple use interface 
to matlab tools for this purpose. I don't  know about mathematica.

However, creating some tools to handle static PDE problems or eigen-value 
problems is not very hard thing. For examples, what's  steps are needed to 
solve eigen-value problem for some PDE in some geometry? We can use, e.g., 
the spectral methods based on Bubnov-Galerkin procedure. 

1)   define the choice of eigen-functions and scale to interval where basis is 
orthogonal. Hints for this come from geometry and from the view of PDE to 
solve.  In Maxima's orthopoly package there is  wide choice for polynomial 
which may be used as a basis.  
2) Creating the spectral matrices need some integration procedure. In simple 
case maxima's functions like integrate, or quad package can be used. It may 
not work in complicated problems.  Instead, we can use Gauss-Legendre 
(LegendreP roots within interval) or Gauss-Raadu (LegendreP extrema + end 
points).  In this case we need collocation points and weights. The examples 
for these routines can be found in attachment. gaulegR(x1,x2,n) for LegendreP 
roots and weights
and gaulegRL(x1,x2,n) for LegendreP extrema + end points.


3)Sometimes we have to reorganize the original polynomial basis to match the 
boundary conditions  . The easiest way is to choose some simple function, say 
F_0, which satisfy the boundary conditions automatically. Then we include F_0 
to the basis and make the Gramm-Schmidt orthogonalization of  the new basis 
set. See further examples in attchment as well.

4)Prepare the matrices for eigen-value problem. Substitute the basis set to 
PDE operator and  integrate 
5)call dgeev on the matrix. 
4 & 5 depends on the task. In fact, the file in attachment is part of bigger 
eigen-value solar dynamo problem. I will present some results about it at the 
meeting in Potsdam "Meridional flow, differential rotation? solar and stellar 
activity", next week. The poster can be found at 
ftp://sftp:sftp at iszf.irk.ru/pipin/postF.pdf 

So my feeling is that creating helpfull tools to handle PDE is not a hard 
problem. Note that it can be extended further  to attack the time-stepping 
problem using ideas from contrib/gentran. I think,  especially for today, 
Maxima has a power! 

Hope it was helpfull! 

rgds
V


-------------- next part --------------
A non-text attachment was scrubbed...
Name: pseudP.mc
Type: text/x-objcsrc
Size: 3735 bytes
Desc: not available
Url : http://www.math.utexas.edu/pipermail/maxima/attachments/20070623/93ba5f51/attachment-0001.bin