Picard’s Method picard.mws

Douglas B. Meade (meade@math.sc.edu)

1 Objectives

• to introduce Picard’s method in a manner accessible to students

• to develop a Maple implementation of Picard’s method, picard

• to use picard to motivate discussion of the existence theory for IVPs

2 New Maple Commands

• for .. do .. od Maple’s primary looping command

• assume command for supplying additional information about a variable or name

• map Maple command to apply a function to all elements of a list, set, or equation (see also map2)

• proc creates a Maple procedure (function) from a sequence of one or more Maple commands

• rsolve Maple command for the solution of a recurrence equation (very similar to dsolve)

• simplify( .., power, symbolic ) forces the symbolic simplification of powers without regard to problems in complex plane

• unapply converts a Maple expression into a Maple function with selected arguments

3 Background

First-order initial value problems of the form \frac{∂}{∂x} y(x) = f(x, y(x)), y(x_{0}) = y_{0} can be rewritten as an integral equation

y(x) = y0 + ∫_{x0}^{x} f(t, y(t)) dt.

This integral formulation can be used to construct a sequence of approximate solutions to the IVP. The basic idea is that given an initial guess of the approximate solution to the IVP, say ϕ_{0}(x), an infinite sequence of functions, {ϕ_{n}(x)} constructed according to the rule

ϕ_{n+1}(x) = y0 + ∫_{x0}^{x} f(t, ϕ_{n}(t)) dt.

That is, the n-th approximation is inserted into the right-hand side of the integral equation in place of the exact solution y(x) and used to compute the ( n + 1)-st element of the sequence. This process is called Picard’s Method.

One of the main applications of Picard’s Method is in the proof of the existence and uniqueness results for first-order initial value problems. This proof will not be considered in this report – hey, we have to leave something for next year’s students to do. This report focuses on gaining an understanding of Picard’s Method through the consideration of a few examples.

4 Example 1

Use Picard’s Method with ϕ_{0}(x) = 1 to obtain the next four successive approximations of the solution to

y’( x)= y(x), y(0) = 1.

4.1 Solution

   > restart;
   > with( plots ):
   > with( DEtools ):

Since no specific choice of phi[0] is given, let’s assume the initial guess is the constant function that satisfies the initial condition:

   > phi[0] := 1;

ϕ_{0} := 1

The next approximation can be computed according to the integral formula. This yields

> phi[1] := 1 + int( phi[0](x), t=0..x );

ϕ_{1} := 1 + x

The same process is used to compute all subsequent approximations. However, in order to keep the same syntax in the integral, it’s necessary to ensure that each ϕ_{n} is a Maple function. The unapply command is used to convert an expression into a function in terms of one or more of the names that appears in the expression. (See the on-line help for unapply for a full description of this command.) Thus, the first approximation is, as a function:

> phi[1] := unapply( phi[1], x );

ϕ_{1} := x → x + 1

The second approximation to the solution can be obtained by repeating this process

> phi[2] := 1 + int( phi[1](t), t=0..x ):
> phi[2] := unapply( phi[2], x );

ϕ_{2} := x → 1 + x + \frac{1}{2} x^{2}

For each of the next two approximate solutions, the two commands will be combined into a single Maple command

> phi[3] := unapply( 1 + int( phi[2](t), t=0..x ), x );

ϕ_{3} := x → 1 + x + \frac{1}{2} x^{2} + \frac{1}{6} x^{3}

> phi[4] := unapply( 1 + int( phi[3](t), t=0..x ), x );

ϕ_{4} := x → 1 + x + \frac{1}{2} x^{2} + \frac{1}{6} x^{3} + \frac{1}{24} x^{4}

A pattern to the sequence of approximate solutions is now obvious. Each iterate adds one new term. The new term for ϕ_{n}(x) is \frac{x^{n}}{n!} . Thus, for each n>=0,

ϕ_{n}(x) = ∑_{k=0}^{n} \frac{x^{k}}{k!} .

These functions are familiar – they are the partial sums of the Maclaurin series for e^{x}.

From our knowledge of power series, we know that these partial sums converge to the exponential function for all real numbers x. This convergence can be demonstrated graphically

> plot( [ exp(x), phi[i](x) $ i=0..4 ], x=0..3,
> title=`Picard Iterates and Limit for y'=y, y(0)=1` );

To conclude this example, we note that the exponential function is a solution to the initial value problem.

> dsolve( { diff( y(x), x ) = y(x), y(0)=1 }, y(x) );

y(x) = e^{x}

5 Example 2

Use Picard’s Method with ϕ_{0}(x) = 0 to obtain the next three successive approximations of the solution to the nonlinear problem y’ = 2 x − y^{2}, y(0) = 0.

5.1 Solution

These approximations could be obtained exactly as in Example 1. However, it will ultimately be more efficient to create a Maple implementation of Picard’s Method. The basic ingredients for Picard’s Method are the right-hand side of the ODE (in standard form), the initial condition, and the initial/previous approximate solution. Here is our version of the method

> picard := proc(f, init, phi)
> local x0, y0;
> x0 := init[1];
> y0 := init[2];
> unapply( y0 + int( f(t,phi(t)), t=x0..x ), x );

> end;

picard := proc(f, init, ϕ)

localx0 , y0 ;

x0 := init_{1} ; y0 := init_{2} ; unapply(y0 + int(f(t, ϕ(t)), t = x0 ..x), x) end

Note that the initial condition, init, is expected to be the list of two elements (i.e., an ordered pair) [x0,y0] and f and phi are assumed to be Maple functions of two and one variable, respectively. The output from picard is the next approximate solution, as a Maple function (so that, e.g., the output is suitable for immediate use in another call to picard.

To use picard to provide the requested information in this example, let

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

f := (x, y) → 2 x − y^{2}

and

> phi[0] := 0;

ϕ_{0} := 0

The first Picard iterate is

> phi[1] := picard( f, [0,0], phi[0] );

ϕ_{1} := (x → x)^{2}

5.1.1 Explanation of Output

This output says that ϕ_{2} is the square of the identity function. Thus, it is simply a fancy way of saying that ϕ_{2}(x) = x^{2}.

The next two iterates are found to be

> phi[2] := picard( f, [0,0], phi[1] );

ϕ_{2} := x → x^{2} − \frac{1}{5} x^{5}

> phi[3] := picard( f, [0,0], phi[2] );

ϕ_{3} := x → x^{2} − \frac{1}{275} x^{11} + \frac{1}{20} x^{8} − \frac{1}{5} x^{5}

The problem simply asks for a plot of these approximate solutions on the interval 0 <= x <=1. The picture is a little nicer if a few more iterates are included. The computation of the next two iterates can be done in a simple for ... do ... od command (see ?do for an explanation of the syntax). For example,

> for n from 0 to 4 do
> phi[n+1] := picard( f, [0,0], phi[n] );
> od;
> n := 'n':

ϕ_{1} := (x → x)^{2}

ϕ_{2} := x → x^{2} − \frac{1}{5} x^{5}

ϕ_{3} := x → x^{2} − \frac{1}{275} x^{11} + \frac{1}{20} x^{8} − \frac{1}{5} x^{5}

ϕ_{4} := x → x^{2}−\frac{1}{5} x^{5}−\frac{1}{1739375} x^{23}+\frac{1}{55000} x^{20}−\frac{87}{374000} x^{17}+\frac{3}{1540} x^{14}−\frac{7}{550} x^{11}+\frac{1}{20} x^{8}

ϕ_{5} := x → x^{2} − \frac{1363623}{3361446550000} x^{29} − \frac{1}{142194993359375} x^{47} + \frac{384503}{123008600000} x^{26}

+ \frac{1}{2104643750000} x^{44} − \frac{8751}{413971250} x^{23} − \frac{7781}{533431525000000} x^{41} + \frac{659}{5236000} x^{20}

+ \frac{194877}{692159930000000} x^{38} − \frac{1769}{2618000} x^{17} − \frac{6915599}{1734042772000000} x^{35} + \frac{1}{308} x^{14}

+ \frac{236231}{5298832000000} x^{32} − \frac{7}{550} x^{11} + \frac{1}{20} x^{8} − \frac{1}{5} x^{5}

The convergence of these polynomials is not as obvious as the first example. Even from the graph it’s difficult to determine the interval on which this integral exists. Including the direction field in the plot shows precisely how these approximate solutions must behave. The following graph provides the essential information for this problem. (The plots on [0, 1] do not show too much difference between the successive approximations. Doubling the interval to [0, 2] shows a little more of the interesting features of these solutions.

> P1 := plot( { phi[i](x) $ i=0..6 }, x=0..2, color=GREEN ):
> P2 := DEplot( diff( y(x),x ) = f(x,y(x)), y, x=0..2, { [0,0]
},
> arrows=LINE, linecolor=BLUE ):
> display( { P1, P2 },
> title=`Picard Iterates and Direction Field:  y'=2x-y^2, y(0)=0`
);

6 Example 3

In previous work it has been demonstrated that \frac{∂}{∂x} y(x) = 3 y(x)^{( \frac{2}{3} )}, y(2) = 0 does not have a unique solution. Show that Picard’s Method, starting with ϕ_{0}(x) = 0 converges to y(x) = 0, whereas Picard’s method beginning with ϕ_{0}(x) = x − 2 converges to the second solution y(x) = (x − 2)^{3}.

6.1 Solution

Each of these cases can be easily computed using picard (defined in Example 2, above). But, the first part is so simple, it’s really a waste of time and energy to use picard or Maple. The claim is that ϕ_{n}(x) = 0 for all x and for all integers n. Thus, these approximate solutions converge to the zero solution.

The successive approximations generated when the initial approximation to the solution is

   > phi[0] := x -> x-2;

ϕ_{0} := x → x − 2

could be computed using picard. The right-hand side of the ODE is

   > f := (x,y) -> 3*y^(2/3);

f := (x, y) → 3 y^{2/3}

and the initial condition is

   > IC := [2,0];

IC := [2, 0]

The next two iterates are

> phi[1] := picard( f, IC, phi[0] );

ϕ_{1} := x → \frac{9}{5} (x − 2)^{5/3}

> phi[2] := picard( f, IC, phi[1] );

ϕ_{2} := x → ∫_{2}^{x} \frac{3}{5} 9^{2/3} 5^{1/3} ((t − 2)^{5/3})^{2/3} dt

6.1.1 Compatibility with Release 3 Release 3 does not have a problem with the rational exponents. This was one of the changes that was made for Release 4. While this seems to be a silly change, there are many situations where this simplification would be inappropriate. The current situation reduces the likelihood that Maple will return an erroneous result.

Note that Maple does not complete the evaluation of the last integral. This is due to the rational exponents (with odd denominators) that appear in the integrand. There are a variety of methods that can be used to force Maple to make the necessary simplifications. The approach adopted here is to make a slight modification to the definition of picard:

> picard2 := proc(f, init, phi)
> local x0, y0;
> x0 := init[1];
> y0 := init[2];
> unapply( y0 + int( simplify( f(t,phi(t)), power, symbolic ),
t=x0..x ), x );
> end:

The first three iterates are

> phi[1] := picard2( f, [2,0], phi[0] );

ϕ_{1} := x → \frac{9}{5} (x − 2)^{5/3}

> phi[2] := picard2( f, [2,0], phi[1] );

ϕ_{2} := x → \frac{27}{95} (x − 2)^{19/9} 9^{2/3} 5^{1/3}

> phi[3] := picard2( f, [2,0], phi[2] );

ϕ_{3} := x → \frac{81}{6175} (x − 2)^{\frac{65}{27}} 27^{2/3} 95^{1/3} 9^{4/9} 5^{2/9}

Observe that the first three approximates all have the form ϕ_{n}(x) = c_{n}(x − 2)^{r_{n}} for constants c_{n} and r_{n}. Explicit formulae for these constants are not apparent from the first three approximate solutions.

An alternate approach to the problem of finding these sequences and, hopefully, lim_{n→∞} ϕ_{n}, is to obtain a recurrence relation between successive terms in each sequence. To determine the appropriate recurrence relations, suppose the n-th iterate is

> phi[n] := x -> c*(x-2)^r;

ϕ_{n} := x → c (x − 2)^{r}

Then, the next iterate will be

> phi[n+1] := picard2( f, [2,0], phi[n] );

ϕ_{n+1} := x → lim_{t→2+} − 3 \frac{(t − 2)^{(2/3 r+1)} c^{2/3}}{\frac{2}{3} r + 1} + 9 \frac{((x − 2)^{r})^{2/3} c^{2/3} (x − 2)}{2 r + 3}

Maple is hesitant about evaluating the limit because r <= 0 has not been excluded. However, examination of the first three iterates together with the above result indicates that the exponents will always be positive.

> assume( rr>0 );

Then, the ( n + 1)-st Picard iterate is

> simplify( subs( r=rr, phi[n+1](x) ) );

9 \frac{((x − 2)^{rr}˜ )^{2/3} c^{2/3} (x − 2)}{2 rr˜ + 3}

This provides a recurrence equation for r_{n}:

> REr := r(n+1) = 1+2/3*r(n);

REr := r(n + 1) = 1 + \frac{2}{3} r(n)

The solution to this recurrence equation, with initial condition r(0) = 1, can be found using rsolve (see ?rsolve)

> SOLNr := rsolve( { REr, r(0)=1 }, { r(n) } );

SOLNr := \{r(n) = −2 (\frac{2}{3})^{n} + 3\}

Thus, it is now clear that

> map( limit, op(SOLNr) , n=infinity );

lim_{n→∞} r(n) = 3

In the same way, the recurrence equation for the coefficients c_{n} is

> REc := c(n+1) = 9/(3+2*r(n))*c(n)^(2/3);

REc := c(n + 1) = 9 \frac{c(n)^{2/3}}{3 + 2 r(n)}

> REc := subs( SOLNr, REc );

REc := c(n + 1) = 9 \frac{c(n)^{2/3}}{9 − 4 (\frac{2}{3})^{n}}

with initial condition c_{0} = 1, is not so easily solved (by hand or by Maple).

> rsolve( { REc, c(0)=1 }, { c(n) } );

rsolve (\{c(0) = 1, c(n + 1) = 9 \frac{c(n)^{2/3}}{9 − 4 (\frac{2}{3})^{n}}\} , \{c(n)\})

It is, nonetheless, still possible to determine the limiting behavior of c_{n}. Note that, since lim_{n→∞} r(n) = 3, lim_{n→∞} \frac{c_{n+1}}{c_{n}^{( \frac{2}{3} )}} = lim_{n→∞} \frac{9}{3+2 r_{n}} = 1. Then,

assuming that c_{n} converges to, say, c_{∞} then c_{∞}^{( \frac{1}{3} )} = 1 and so c_{∞} = 1, since these coefficients are real-valued.

This completes the proof that this sequence of Picard iterates converges, with a limit of y(t) = lim_{n→∞} ϕ_{n}(t) = (t − 2)^{3}.

6.1.2 Point to Ponder

• Verify that y(t) = (t − 2)^{3} is a solution to the IVP.

The result that is obtained from these computations is that this IVP does not have a unique solution. This conclusion is consistent with the existence and

uniqueness theory (since \frac{∂}{∂y} f(x, y) = \frac{2 y^{(− \frac{1}{3} )}}{3} is not defined – not to mention continuous – at y = 0.)

7 Example 4

Use Picard’s method with ϕ_{0}(x) = 1 to compute the next six successive approximations of the solution to the nonlinear problem y’ = 1 + y^{2}, y(0) = \frac{1}{2}. Create plots of the each Picard iterate superimposed on the direction field. What can you say about the existence of a solution to this IVP?

7.1 Solution

The data required by picard is given as in the previous examples:

> f := (x,y) -> 1+y^2;
> IC := [ 0, 1/2 ];
> phi[0] := 1;

f := (x, y) → 1 + y^{2}

IC := [0, \frac{1}{2}]

ϕ_{0} := 1

The computation of the Picard iterates is straightforward. Only the first four iterates are displayed, for the sake of saving a few pages of output.

> for n from 0 to 5 do
> phi[n+1] := picard( f, IC, phi[n] );
> od:
> phi[1](t);

\frac{1}{2} + 2 t

> phi[2](t);

\frac{1}{2} + \frac{5}{4} t + t^{2} + \frac{4}{3} t^{3}

> phi[3](t);

\frac{1}{2} + \frac{5}{4} t + \frac{16}{63} t^{7} + \frac{4}{9} t^{6} + \frac{13}{15} t^{5} + \frac{23}{24} t^{4} + \frac{41}{48} t^{3} + \frac{5}{8} t^{2}

> phi[4](t);

\frac{1}{2} + \frac{5}{4} t + \frac{24121}{60480} t^{9} + \frac{256}{59535} t^{15} + \frac{5491}{10752} t^{8} + \frac{64}{3969} t^{14} + \frac{10457}{16128} t^{7} + \frac{1808}{36855} t^{13} + \frac{4157}{5760} t^{6}

+ \frac{11}{105} t^{12} + \frac{223}{320} t^{5} + \frac{3208}{17325} t^{11} + \frac{29}{48} t^{4} + \frac{10349}{37800} t^{10} + \frac{35}{48} t^{3} + \frac{5}{8} t^{2}

The convergence of the constant, linear, and quadratic terms are apparent from these results. One additional term seems to be correct after each Picard iterate, but the total number of terms doubles with each iteration.

> P1 := DEplot( diff( y(x), x ) = f(x,y(x)), y, x=0..2, {IC}, y=0..10,
> arrows=LINE, linecolor=BLUE ):
> P2 := plot( { phi[i](x) $ i=0..6 }, x=0..2, y=0..10, color=GREEN
):
> display( [ P1, P2 ],
> title=`Direction Field and Picard Iterates:  y'=1+y^2, y(0)=1/2`
);

7.1.1 Point to Ponder Note how close the last few Picard iterates are to the exact solution.

• Compare this plot with the one produced when the order of the plots is reversed in the display command.

Animating the iterates clearly demonstrates the convergence of the Picard iterates to the solution of this IVP (at least on some interval).

> OPTIONS := y=0..10, color=GREEN :
> for i from 0 to 6 do
> p.i := display( [ P1, plot( phi[i](x), x=0..2, OPTIONS ) ] );
> od:

> i := 'i':
> display( [ seq( p.i, i=0..6 ) ], insequence=true, view=[0..2,0..10],
> title=`Direction Field and Picard Iterate for y' = 1 + y^2, y(0)=1/2`
);

7.1.2 Suggestions for “Plotting an Animation” (optional) Note that the above animation is not easily printed. The new feature of Release 4 that displays an array of plots is the best option at the present time (see ?plots,display for more information on this feature).

> display( array( 1..2, 1..3, [ [p1, p2, p3], [p4, p5, p6] ] )
);

Note that some options are ignored and/or overwritten when using display. This is particularly noticeable, as in the previous example, where an array of plots is displayed. In this case it is the view option that cannot be specified. (It does no good to include view=[0..2,0..10] in OPTIONS or the final display command (this option is used appropriately in the creation of the animation).

The graphical information also suggests that the solution may not exist for all x>0. The solution seems to have a vertical asymptote, x = x_{crit}, for some x_{crit} slightly larger than 1. In fact, this problem is not difficult to solve by hand. The exact solution is y(x) = tan(x + arctan(\frac{1}{2})) which has a vertical asymptote at x = x_{crit} where the argument to the tangent function is \frac{π}{2} ; that is:

   > x[crit] := solve( x+arctan(1/2) = Pi/2, x );

x_{crit} := −arctan(\frac{1}{2}) + \frac{1}{2} π

   > evalf( x[crit] );

1.107148718

This is consistent with the graphical evidence that has been presented.

8 Conclusion

In this project we have developed a Maple implementation of Picard’s method and illustrated the use of this procedure for four different initial value problems. These examples have demonstrated the convergence of Picard iterates to solutions of the IVP. The third example is most interesting in that even though the individual successive approximations could not be determined explicitly, the limit of the sequence could still be resolved. (This example was also interesting in that it was an example of an IVP with multiple solutions.) The fourth example is an additional problem that is selected to illustrate the convergence of the Picard iterates when the solution exists only on a finite interval.

9 Copyright

picard.mws – a Maple worksheet for Picard’s method (for use with Chapter 1, Project B, in Nagle and Saff) Prepared by: Douglas B. Meade (7 March 1996) as part of the Maple Supplement for Nagle/Saff (Addison-Wesley) Updated to Maple V, Release 4, by Douglas B. Meade (27 December, 1996) For more information, please contact the author at one of following addresses: Department of Mathematics, Univ. of South Carolina, Columbia, SC 29208 (803) 777-6183 meade@math.sc.edu (copyright, Douglas B. Meade, 1996)

15