Lab 12 -- Picard Iteration

---

Kris's Lab 12 Walkthrough | Lab 12 handout | My Picard Method Diagram

We're here going to run through Picard iteration.  Picard iteration is a special kind of fixed point iteration.  We call x a fixed point of of a function F, iff. x = F(x).  Suppose you define a sequence by 
  x_n+1 = F(x_n), 
  x₁ = [some guess at the fixed point].
Often you will find that x_n converges to a fixed point of F. The process of taking the successive terms of such a sequence is called iteration. We're going to apply this iterative idea to differential equations, and we come up with the Picard method. Basically, we're going to apply fixed point iteration to a whole differential equation.

How do we do that? To start, have a look at this diagram on the Picard method. It might not make sense yet, but don't worry. It might help to keep a printout of this at hand as you go through this walkthrough; I'm going to refer to notation from that diagram in what follows. So let's use this notation:

• We have a differential equation, call it [*]
    y'(t) = G(t, y(t));   
    y(0) = y₀.
  Our goal is to use the Picard method to try to find a solution to it. How does the Picard method work?
• Well, first notice that 
    y(t) - y₀  =  I(t, y(t)), 
  where 
    I(t, y(t)) =  the definite integral of G(s, y(s)) with respect to s, from s = 0 to s = t.  
  (This is just because y is an antiderivative of y'.) 
• In the Picard method, we rewrite this statement as our fixed point equation, y(t) = y₀ + I(t, y(t)). Let F(t, y(t)) denote the right-hand side of the new fixed point equation. So F(t, y(t)) = y₀ + I(t, y(t)), and y(t) = F(t, y(t)).
• Note that any solution to our fixed point equation just stated would also solve our original differential equation, [*]. Remember that fixed point iterations often converge to a solution of the fixed point equation.
• Thus, our fixed point sequence will be defined by y_k+1(t) = F(t, y_k(t)). 
• Finally, let "starting condition" for this sequence be  y₁(t) = y0.

The Picard method basically works by using y₁(t) to get y₂(t), and using y₂(t) to get y₃(t), and so on. The hope is that y_k(t) will get closer to a solution to the fixed point equation (and thus to a solution of our differential equation, [*]) as k gets bigger. 

We're next going to actually perform these fixed point iterations in MATLAB for a sample equation, creatively enlisting the help of ODE23. You might want to glance back at the Picard method diagram now. Everything there should make sense to you now, except for where and why ODE23 fits into this whole process.

How you can do this with Matlab, and why you can do it this way: A specific example

In Section 2, #1a, you're asked to compute the Picard iterate of a given differential equation, equation (3):
  y' = G(t, y(t)) = t² + y²
  y(0) = y0 = seed / 10 (My section: Use seed = 18. Others: substitute your own seed)
We want the Picard iterate computed for t in [0, 4].

Now, remember that y(t) = F(t, y(t)) = y₀ + I(t, y(t)), where I(t, y) [defined as it is above in this document] is the integral you see in equation (6) in the Lab 12 handout. Also, as always, the Picard iterates are defined by y_k+1(t) = F(t, y_k(t)) and y₁(t) = y₀ for any t.

So let's start: Let's compute (as they ask us to do in #1a) y₂(t). Using our formulas, we get:

y₂(t) = F(t, y₁(t)) = y₀ + I(t, y₁(t)) = seed / 10 + I(t, y(t)) = 1.8 + I(t, y(t)).

But how do we find I(t, y(t)) (the integral of G(t, y(t)) from 0 to t)? Well, the solution set to the integral is the same as the solution set to the differential equation that comes from it:
  dI / dt = G(t, y(t)) = t^2 + y^2.
  I(0, y) = 0  (any integral from 0 to 0 is 0).
This is a differential equation. Thus, we can approximate the solution set, I(t, y) with ode23. And that will give us the integral.

Now, how to go about the mechanics of that? ODE23 needs a derivative function, and an initial condition for the solution (i.e., the value of I(0). The problem is that this differential equation in I(t, y) has a derivative in terms of t and y, not in terms of t and I. But the ODE23 syntax is restrictive,  and requires that the parameters passed to the derivative function are only the values of t and I. So we now have to find a way for our derivative function to know the values of t and y without being passed these values as parameters. That's why we use the "global" command -- we set up y and t as global variables that can be accessed from inside any function.