Picard Iteration for solving ODEs/PDEs

The slides from the Maple Summer Workshop 2004 talk are here.

The Maple worksheet for the Springy-Pendulum is here: Spring-Pendulum-AMP2.mws.

The files for the burger example are here:

Input for C++ Code	burger.in
Output for 100 steps at h=0.1	burger.out
Maple worksheet to load the data and display it	burger.mws

This is the animation of the solar system using the Modified Picard method.

Finally, the code to run the C++ examples is pdepoly.tar.gz. You will need MPI to use this. See below to use this.

I have developed a parallel algorithm that will produce Maclaurin Polynomial solutions to any initial value ODE or PDE using Picard iteration as discussed in the paper "Implementing the Picard Iteration" appearing in Neural, Parallel & Scientific Computations (March 1996 pp. 97-112) by G. Edgar Parker and James Sochacki. The basic idea is to convert the ODE or PDE to polynomial form using elementary substitutions and then solving the resulting system of ODE's or PDE's using Picard iteration. My code outputs the polynomials in a form that can be read into Maple and then one can plot the polynomials to visualize the solution. I give an example below on how to use the code, but you may want to read my paper "A General ODE and PDE Solver Using Picard's Method" for a more detailed description. Please email me any comments on your attempts at using this software.

The code requires some version of MPI to run. I no longer maintain a version the runs without the use of MPI. I recommend using MPICH or LAM, as I have tested my code on both systems and it works with them.

Input

The input is in the form of :


	sp[<number of space variables>];

	var[<number of equations>];

	var[1]':=<equation>;

	var[1](0):=<initial value>;

	...

	var[<number of equations>]':=<equation>;

	var[<number of equations>](0):=<initial value>;

We will use t as the variable in ODE's and as the initial value variable in PDE's. The other variables in PDE's will be referred to as space variables. The t variable will also be referred to as time. We will denote a t derivative in both ODE's and PDE's as '. The ODE or PDE must be independent of t. For example:


	U'=Ux

In this example, x is a space variable and U is an equation variable.

Now, an equation can consist of equation vaiables raised to any power, as well as any space variable raised to a power. Both can have an associated coefficient in the front, and you can just have a number such as 1.0.
The initial value is the same as an equation except that it can't have equation variables in its definition.
As an example, consider:

Y' = cos(Y) + sin(t); Y(0) = 0

Remember we have to transform this equation to polynomial form in all variables and the equation has to be independent of t (autonomous). Making the substitution U=cos(Y),V=sin(t),W=sin(Y) and Z=cos(t) gives:


U'= -Z(U+V) ; U(0) = 1

V'= Z ; V(0) = 0

W'= U*(U+V) ; W(0) = 0

Z'= -V ; Z(0) = 1

The format for this system in my code is:


	sp[0];

	var[5];

	var[1]':=var[2]+var[3];

	var[1](0):=0;

	var[2]':=-1.0var[4]*var[2]+-1.0var[4]*var[3];

	var[2](0):=1;

	var[3]':=var[5];

	var[3](0):=0;

	var[4]':=var[2]*(var[2]+var[3]);

	var[4](0):=0;

	var[5]':=-1.0var[3];

	var[5](0):=1;

Also, the equation entry allows for several additions:

An equation inside a () like: var[1]*(var[2]+var[3]).

The partial derivative using D(): D(var[1]*sp[2]^3,sp[1]), which takes the derivative of U*y^3 with respect to x. You can also take the derivative of a derivative and so on.

Output

The output returned to you cantains any errors plus the output used in Maple V Release 5 and later versions. It is in the form of :


	p[<equation number>,<time step, 0 is initial value>]:=<answer>;

	...

	up to number of equations*number of time steps.

You can load this in Maple using the read command:


	read("data.mpl");

Date: July 12, 2004.