Picard's Iteration Method

What is the Purpose of Picard’s Iteration Method?

The real power of Picard’s Iteration Method is that it can be used to prove the existence and uniqueness theorem for nonlinear initial value problems. To do this, you have to determine under what conditions the successive approximations converge. Here is the theorem (from Coddington, An Introduction to Ordinary Differential Equations)

Let f be a continuous real-valued function in the rectangle 
R:|x-(Subscript[x, o])|≤ a, |y-(Subscript[y, 0])|≤b, a,b>0
and let
|f(x,y)|≤ M
for all (x,y)in R. 

Further suppose f satisfies a Lipshitz condition 
|f(x,(Subscript[y, 1]))-f(x,(Subscript[y, 2]))|≤ K|(Subscript[y, 1])-(Subscript[y, 2])|
with constant K in region R. 

Then the successive approximations
(Subscript[∅, 0])(x)=(Subscript[y, 0]), (Subscript[∅, k+1])(x)=(∫_(Subscript[x, 0])^x)f(s,(Subscript[∅, k])(s))ds
converge on the the interval 
I:|x-(Subscript[x, o])|≤ min{a,b/M}
to a solution of the initial value problem
y'=f(x,y),y((Subscript[x, 0]))=(Subscript[y, 0])
on I

Notice this is slightly different than the existence and uniqueness theorem in our text.

Example 1

Clear[f,ϕ,y0,t0]

t0=0;
y0=1;
f[t_,y_]=y+1-t;

ϕ[0,t_]=y0;
ϕ[n_,t_]:=ϕ[n,t]=y0+(∫_t0^t)f[s,ϕ[n-1,s]]\[DifferentialD]s

sequence=Table[ϕ[n,t],{n,0,12}];
TableForm[sequence]

Plot[sequence,{t,-10,10},PlotRange->{All,{-10,20}},BaseStyle->Large,AxesLabel->{"t","y"}]

Example 2

Clear[f,ϕ,y0,t0]

t0=0;
y0=1;
f[t_,y_]=-y-t;

ϕ[0,t_]=y0;
ϕ[n_,t_]:=ϕ[n,t]=y0+(∫_t0^t)f[s,ϕ[n-1,s]]\[DifferentialD]s

sequence=Table[ϕ[n,t],{n,0,12}];
TableForm[sequence]

Plot[sequence,{t,0,10},PlotRange->{All,{-5,5}},BaseStyle->Large,AxesLabel->{"t","y"}]

Example 3

Clear[f,ϕ,y0,t0]

t0=0;
y0=1;
f[t_,y_]=-y+Cos[t/3];

ϕ[0,t_]=y0;
ϕ[n_,t_]:=ϕ[n,t]=y0+(∫_t0^t)f[s,ϕ[n-1,s]]\[DifferentialD]s

sequence=Table[ϕ[n,t],{n,0,12}];
TableForm[sequence]

Plot[sequence,{t,0,10},PlotRange->{All,{-5,5}},BaseStyle->Large,AxesLabel->{"t","y"}]

Example 4 Nonlinear

For nonlinear, the integrals are generally difficult to evaluate, and the integrand need to be expanded in Taylor series about (Subscript[t, 0]). This method is chosen so we can still get a function, abeit an approximation to the function. You could of course use other approximations such as Fourier series if you like.

This is a terrible way to numerically solve an IVP.

Due to the computational effort required, we will only store the previous result and not all the (Subscript[ϕ, k])(t).

Clear[f,ϕ,y0,t0]

t0=-0.5;
y0=1.0;
f[t_,y_]=y-(t^2)+Sin[t y]t;

ϕ[t_]=y0;

max=100;
k=1;
While[k≤max,
integrand=Normal[Series[f[s,ϕ[s]],{s,t0,70}]];
ϕ[t_]=y0+(∫_t0^t)integrand\[DifferentialD]s;
k=k+1;
]

pPicard=Plot[ϕ[t],{t,-1,5},PlotRange->{All,{0,11}},BaseStyle->Large,AxesLabel->{"t","y"},PlotStyle->Thick]

Here we see Taylor’s Method of Order 2 is much better.

Clear[f,t,y]
f[t_,y_]=y-(t^2)+Sin[t y]t;

h=0.1;
t[0]=-0.5;
y[0]=1.0;
t[i_]:=t[i]=t[i-1]+h
y[i_]:=y[i]=y[i-1]+h*f[t[i-1],y[i-1]]+((h^2)/2)(Derivative[1,0][f][t[i-1],y[i-1]]+Derivative[0,1][f][t[i-1],y[i-1]]*f[t[i-1],y[i-1]])

taylor2=Table[{t[k],y[k]},{k,0,50}];
taylor2plot=ListPlot[taylor2,PlotStyle->{Blue},AxesLabel->{"t","y"}];

Show[pPicard,taylor2plot]