Maple_A2.mw

Appendix A2. Picard Iteration

The following first order initial value problem

y'(t) = f ( t , y(t) )  ,  y(0) = 0

can be reformulated as the an integral equation

y(t) = int(f(s, y(s)), s = 0 .. t) .

See Ledder, Appendix A.2.

Assume that f (t, y) and its y partial are continuous on a rectangle containing the point (0,0). Then starting with an initial "guess" for a solution

y[0](t) = g(t) ,

the recursive relation

y[n](t) = int(f(s, y[n-1](s)), s = 0 .. t)

generates a sequence of functions converging to the solution. These functions are called the Picard iterates generated by g(t).

Maple can generate some of these iterates, use a for..do loop.

The following example is taken from Ledder, see Chapter A.2, Exercise 2.

Example  Obtain the first three Picard iterates for the following IVP, generated by the function g(t) = 0.

diff(y(t), t) = t+t^2*y(t)  ,  y(0) = 0

Sketch their graphs and the graph of the exact solution.

Start with the definition of the function f.

> f := (t,y) -> t + t^2*y;

f := proc (t, y) options operator, arrow; t+t^2*y end proc

Use it to define the diffenential equation.

> DE := diff(y(t),t) = f(t,y(t));

DE := diff(y(t), t) = t+t^2*y(t)

Compute the first three iterates.

> Y[0] := t -> 0:
for n from 1 to 3

 do

   int(f(s,Y[n-1](s)),s=0..t):

   Y[n] := unapply(%,t):

 end do:

unassign('n');

The formulas for the initial guess and the three iterates are displayed below.

> for n from 0 to 3 do y[n](t) = Y[n](t) end do; unassign('n');

y[0](t) = 0

y[1](t) = 1/2*t^2

y[2](t) = 1/2*t^2+1/10*t^5

y[3](t) = 1/80*t^8+1/10*t^5+1/2*t^2

The actual solution formula is obtained next. It is not pretty.

> soln := dsolve( {DE, y(0)=0} );

soln := y(t) = 1/10*exp(1/6*t^3)*9^(2/3)*(t^3*WhittakerM(1/3, 5/6, 1/3*t^3)+5*WhittakerM(4/3, 5/6, 1/3*t^3))/(t*(t^3)^(1/3))

The plot of the solution, red and thick, and the three iterates, blue, green, black:

> plot( [rhs(soln),Y[n](t)$n=1..3], t=0..2, y=0..5,
      color=[red,blue,green,black], thickness=[2,1$3]);

[Plot]

The next output reveals that the Picard iterates are the partial sums to the Taylor series for the solution.

> dsolve( {DE, y(0)=0}, y(t), type=series, order=9);

y(t) = (series(1/2*t^2+1/10*t^5+1/80*t^8+O(t^9),t,9))

Ask the FunctionAdvisor about the Whittaker functions.

> FunctionAdvisor(Whittaker);

The 2 functions in the "Whittaker" class are:

[WhittakerM, WhittakerW]

This is enough to get to a Help page. Go see what they are all about.