Preface
This section expands Picard's iteration process to systems of ordinary differential equations in normal form (when the derivative is isolated). Picard's iterations for a single differential equation \( {\text d}x/{\text d}t = f(t,x) \) was considered in detail in the first tutorial (see section for reference). Therefore, our main interest would be to apply Picard's iteration to systems of first order ordinary differential equations in normal form (which means that the highest derivative is isolated)
\[
\begin{cases}
\dot{x}_1 \equiv {\text d}x_1 /{\text d}t &= f_1 (t, x_1 , x_2 , \ldots , x_n ) , \\
\dot{x}_2 \equiv {\text d}x_2 /{\text d}t &= f_2 (t, x_1 , x_2 , \ldots , x_n ) , \\
\quad \vdots & \quad \vdots \\
\dot{x}_n \equiv {\text d}x_n /{\text d}t &= f_n (t, x_1 , x_2 , \ldots , x_n ) .
\end{cases}
\]
This explicit notation cannot be called compact and informative. Hence, we introduce n-dimensional vectors of unknown variables and slope functions in order to rewrite the above system in compact form. So
the system of differential equations can be written in vector form
\[
\frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}(t, {\bf x} ) ,
\]
where
\[
{\bf x} (t) = \begin{bmatrix} x_1 (t) \\ x_2 (t) \\ \vdots \\ x_n (t) \end{bmatrix} ,
\qquad {\bf f} (t, x_1 , x_2 , \ldots , x_n ) = \begin{bmatrix} f_1 (t, x_1 , x_2 , \ldots , x_n ) \\ f_2 (t, x_1 , x_2 , \ldots , x_2 ) \\ \vdots \\
f_n (t, x_1 , x_2 , \ldots , x_n ) \end{bmatrix}
\]
are n-column vectors. Note that the above system of equations contains
n dependent variables \( x_1 (t), x_2 (t) , \ldots ,
x_n (t) \) while the independent variable is denoted by t,
which may be associated with time. In engineering and physics, it is a custom
to follow Isaac Newton and denote a derivative with respect to time variable t by dot:
\( \dot{\bf x} = {\text d}{\bf x} / {\text d} t. \)
When input function f satisfies some general conditions (usually
required to be Lipschitz
continuity in the dependent variable \( {\bf x} \) ), then
Picard's iteration converges to a solution of the Volterra integral equation
in some small neighborhood of the initial point. Thus, Picard's iteration is
an essential part in proving existence of solutions for the initial value
problems.
Note that Picard's iteration is not suitable for numerical calculations. The reason is not only in slow convergence, but mostly it is impossible, in general, to perform explicit integration to obtain next iteration. Although in case of polynomial input function, integration can be performed explicitly (especially, with the aid of a computer algebra system), the number of terms quickly grows as a snow ball. While we know that the resulting series converges eventually to the true solution, its range of convergence is too small to keep many iterations. In another section, we show how to bypass this obstacle.
If an initial position of the vector \( {\bf x} (t) \) is known, we get an initial value problem:
\begin{equation} \label{EqPicard.1}
\frac{{\text d} {\bf x}}{{\text d}t} = {\bf f}(t, {\bf x} ) , \qquad
{\bf x} (t_0 ) = {\bf x}_0 ,
\end{equation}
where \( {\bf x}_0 \) is a given column vector.
Many brilliant mathematicians participated in proving the existence of a
solution to the given initial value problem more than 100 years ago. Their
proof was based on what is now called the Picard's iteration, named after the
French mathematician Charles Émile Picard (1856--1941) whose theories did much to advance
research in analysis, algebraic geometry, and mechanics.
The Picard procedure is actually a practical extension of the Banach fixed point theorem, which is applicable to continuous contractive function. Since any differential equation involves an unbounded derivative operator, the fixed point theorem is not suitable for it. To bypass this obstacle, Picard suggested to apply the (bounded) inverse operator L-1 to the derivative one \( \texttt{D} . \) Recall that the inverse \( \texttt{D}^{-1} , \) called in mathematical literature as «antiderivative», is not an operator because it assigns to every input-function infinite many outputs. To restrict its output to a single one, we consider the differential operator on the set of functions (which becomes a vector space only when the differential equation and the initial condition are all homogeneous) with a specified initial condition f(x0) = y0. So the derivative operator on this set of functions we denote by L, and its inverse is a bounded integral operator. The first step in deriving Picard's iterations is to rewrite the given initial value problem in equivalent (this is true when the slope function f is continuous in Lipschitz sence) form as Volterra integral equation of second kind:
\begin{equation} \label{EqPicard.2}
{\bf x} ( t ) ={\bf x}_0 + \int_{t_0}^t {\bf f} (s, {\bf x} (s) ) \,{\text d}s .
\end{equation}
This integral equation is obtained upon integration of both sides of the
differential equation \(
{\text d}{\bf x} / {\text d} t = {\bf f}(t, {\bf x}), \) subject to the initial condition.
The equivalence follows from the Fundamental Theorem of Calculus. It suffices to find a continuous function \( {\bf x}(t) \) that satisfies the integral equation within the interval \( t_0 -h < t < t_0 +h , \) for some small value \( h \) since the right hand-side (integral) will be continuously differentiable in \( t . \)
Now we apply a technique known as Picard iteration to construct the required solution:
\begin{equation} \label{EqPicard.3}
{\bf x}_{m+1} ( t ) ={\bf x}_0 + \int_{t_0}^t {\bf f} (s, {\bf x}_m (s) ) \,{\text d}s , \qquad m=0,1,2,\ldots .
\end{equation}
The initial approximation is chosen to be the initial value
(constant): \( {\bf x}_0 (t) \equiv {\bf x}_0
.\) (The sign ≡ indicates the values are identically
equivalent, so this function is the constant).
When input function f(t, x) satisfies some (sufficient)
conditions, it can be shown that this iterative sequence converges to
the true solution uniformly
on the interval \( [t_0 -h, t_0 +h ] \) when
h is small enough.
Actually, the necessary and sufficient conditions for existence of a solution
to a vector initial value problem are still waiting for their discovery.
Of course, we can differentiate the recursive integral relation to obtain a sequence of the initial value problems
\begin{equation} \label{EqPicard.4}
\frac{{\text d}{\bf x}_{m+1} ( t )}{{\text d}t} = {\bf f} (t, {\bf x}_m (t) , \qquad {\bf x}_m (t_0 ) = {\bf x}_0 , \qquad m=0,1,2,\ldots .
\end{equation}
that is equivalent to the original recurrence \eqref{EqPicard.3}. However,this recurrence relation \eqref{EqPicard.4} involves the derivative operator
\( \texttt{D} = {\text d}/{\text d}t . \) Since \( \texttt{D} \) is an unbounded operator, it is not possible to prove convergence of iteration procedure \eqref{EqPicard.4}. On the other hand, the convergence of \eqref{EqPicard.3} involving a bounded integral opeprator can be accomplished directly (see Tutorial I).
Example 1:
■
Example 2:
Consider the initial value problem
Although, the next iteration Mathematica can handle by invoking special functions,
\[
\frac{{\text d}y}{{\text d}x} = \sqrt{\frac{1}{4} + 2\,y} , \qquad y(0) = 1.
\]
We apply Picard's iteration directly by solving the fixed point problem
\[
\phi_{n+1} (x) = 1 + \int_0^x \sqrt{\frac{1}{4} + 2\,\phi_n
(t)}\,{\text d}t , \qquad n=0,1,2,\ldots .
\]
If we start with the initial condition ϕ0 = 1, we almost immediately come to a problem of explicit integration. Indeed, we calculate first the few terms:
\begin{align*}
\phi_0 &= 1 ,
\\
\phi_1 (x) &= 1 +\int_0^x \sqrt{\frac{1}{4} + 2}\,{\text d}t = 1 + \frac{3}{2}\, x ,
\\
\phi_2 (x) &= 1 +\int_0^x \sqrt{\frac{1}{4} + 2\,\phi_1 (t)}\,{\text d}t = 1+ \int_0^x \sqrt{\frac{1}{4} + 2 +3\,t}\,{\text d}t = \frac{x}{\sqrt{3}}\, \sqrt{3 + 4\,x} + \frac{1}{4} \left( 1 + \sqrt{9 + 12\,x} \right) ,
\end{align*}
because Sage helps:
Simplify[1 +
Assuming[x > 0,
Integrate[
Sqrt[t Sqrt[1 + (4 t)/3] + 1/4 (1 + Sqrt[9 + 12 t])], {t, 0, x}]]
it will not overcome the further step.
To apply Picard's method, we convert the given initial value problem to the system of differential equations upon introducing auxiliary dependent variables:
\begin{align*}
y_1 (x) &= y(x),
\\
y_2 (x) &= \sqrt{\frac{1}{4} + 2\,y(x)} ,
\\
y'_1 &= y_2 (x) ,
\\
y'_2 (x) &= \frac{1}{2\,\sqrt{\frac{1}{4} + 2\,y(x)}} \, 2\, y' (x) = 1 .
\end{align*}
We get the vector equation with polynomial right hand side:
\[
\frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (x) \\ y_2 (x)
\end{bmatrix} = \begin{bmatrix} y_2 (x) \\ 1
\end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0)
\end{bmatrix} = \begin{bmatrix} 1 \\ \frac{3}{2}
\end{bmatrix} .
\]
Integrating the vector equation above allows us to rewrite the given initial value problem in
equivalent vector form:
\[
\begin{bmatrix} y_1 (x) \\ y_2 (x) \end{bmatrix} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} y_2 (s) \\ 1 \end{bmatrix} {\text d} s .
\]
The Picard's iteration scheme will be
\[
\begin{bmatrix} \phi_1 (x) \\ \phi_2 (x) \end{bmatrix}_{(m+1)} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} \phi_2 (s) \\ 1 \end{bmatrix}_{(m)} {\text d} s , \qquad m=0,1,2, \ldots .
\]
Since Picard's iteration always starts with the initial condition, we get
\[
\begin{bmatrix} \phi_1 (x) \\ \phi_2 (x) \end{bmatrix}_{(1)} =
\begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x
\begin{bmatrix} \frac{3}{2} \\ 1 \end{bmatrix} {\text d} s =
\begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \begin{bmatrix}
\frac{3}{2}\, x \\ x \end{bmatrix} = \begin{bmatrix}
1+\frac{3}{2}\, x \\ \frac{3}{2} + x \end{bmatrix} .
\]
Next term will be
\[
\begin{bmatrix} \phi_1 (x) \\ \phi_2 (x) \end{bmatrix}_{(2)} = \begin{bmatrix} 1 \\ \frac{3}{2} \end{bmatrix} + \int_0^x \begin{bmatrix} \frac{3}{2} +s \\ 1 \end{bmatrix} {\text d} s = \begin{bmatrix} 1+ \frac{3}{2}\, x + \frac{x^2}{2} \\ \frac{3}{2} +x \end{bmatrix} .
\]
The second component, ϕ2, remains the same with
every iteration. This component is used to generate the first
component, ϕ1, as well. As such, every following iteration will be the
same, and the true solution becomes
\[
y(x) = 1 +\frac{3}{2}\, x + \frac{x^2}{2} .
\]
Its derivative is
\[
y' (x) = \frac{3}{2} + x \qquad \Longrightarrow \qquad \frac{1}{4} + 2\, y(x) = \left( \frac{3}{2} + x \right)^2 .
\]
■
Example 3:
Simple pendulum equation
Consider the initial value problem for simple pendulum equation
\[
\frac{{\text d}^2 \theta}{{\text d} t^2} + \sin \theta = 0 , \qquad \theta (0) = \frac{\pi}{4} , \quad \dot{\theta} (0) = \frac{1}{10} ,
\]
where dot stands for the derivative with respect to time
variable t.
We integrate the pendulum equation above twice, and use initial condition
provided. We may try to apply Picard's iteration scheme
\[
\phi_{m+1} (t) = \frac{\pi}{4} + \frac{1}{10}\, t - \int_0^t \left( t-s \right) \sin \phi_m (s) \,{\text d}s , \qquad m= 1,2, \ldots .
\]
Although we can find the second term
\[
\phi_{2} (t) = \frac{\pi}{4} + \frac{1}{10}\, t - \int_0^t \left( t-s \right) \sin \left( \frac{\pi}{4} + \frac{1}{10}\,s \right) {\text d}s = \frac{\pi}{4} + \frac{1}{10}\, t + 5\sqrt{2} \left( 10 +t \right) - 100\,\sin \frac{5\pi + 2t}{20} .
\]
Simplify[Pi/4 + t/10 +
Assuming[t > 0, Integrate[(t - s)*Sin[Pi/4 + s/10], {s, 0, t}]]]
next term ϕ3(t) is impossible to determine, even for such powerful computer algebra system as Mathematica. We can bypass this obstacle by introducing auxiliary variables:
\begin{align*}
y_1 (t) &= \theta (t),
\\
y_2 (t) &= \dot{y}_1 = \dot{\theta} ,
\\
y_3 (t) &= \sin \theta (t) ,
\\
y_4 (t) &= \cos \theta (t) .
\end{align*}
Then the pendulum equation and corresponding initial conditions will be
\[
\frac{\text d}{{\text d}t} \begin{bmatrix} y_1 (t) \\ y_2 (t) \\ y_3 (t) \\ y_4 (t) \end{bmatrix} = \begin{bmatrix} y_2 (t) \\ -y_3 (t) \\ - y_2 (t)\,y_4 (t) \\ y_2 (t)\,y_3 (t) \end{bmatrix} , \qquad \begin{bmatrix} y_1 (0) \\ y_2 (0) \\ y_3 (0) \\ y_4 (0) \end{bmatrix} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix}
\]
Since the above system of four differential equations has polynomial slope part, we expect to perform all steps of Picard's iteration scheme:
\[
\begin{bmatrix} \phi_1 (t) \\ \phi_2 (t) \\ \phi_3 (t) \\ \phi_4 (t) \end{bmatrix}_{(m+1)} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \int_0^t
\begin{bmatrix} \phi_2 (t) \\ -\phi_3 (t) \\ - \phi_2 (t)\,\phi_4 (t) \\ \phi_2 (t)\,\phi_3 (t) \end{bmatrix}_{(m)} {\text d} s , \qquad m=0,1,2,\ldots .
\]
The first term is
\[
\begin{bmatrix} \phi_1 (t) \\ \phi_2 (t) \\ \phi_3 (t) \\ \phi_4 (t) \end{bmatrix}_{(1)} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \int_0^t
\begin{bmatrix} \frac{1}{10} \\ -\frac{1}{\sqrt{2}} \\ - \frac{1}{10} \,\frac{1}{\sqrt{2}} \\ \frac{1}{10} \,\frac{1}{\sqrt{2}} \end{bmatrix} {\text d} s = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + t \begin{bmatrix} \frac{1}{10} \\ -\frac{1}{\sqrt{2}} \\ - \frac{1}{10} \,\frac{1}{\sqrt{2}} \\ \frac{1}{10} \,\frac{1}{\sqrt{2}} \end{bmatrix} .
\]
Next,
\[
\begin{bmatrix} \phi_1 (t) \\ \phi_2 (t) \\ \phi_3 (t) \\ \phi_4 (t) \end{bmatrix}_{(2)} = \begin{bmatrix} \frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} + \int_0^t
\begin{bmatrix} \frac{1}{10} - \frac{s}{\sqrt{2}}
\\
-\frac{1}{\sqrt{2}} + \frac{s}{10\,\sqrt{2}}
\\
- \left( \frac{-s}{\sqrt{2}} + \frac{1}{10} \right)
\left( \frac{1}{\sqrt{2}} + \frac{s}{10\,\sqrt{2}} \right)
\\
\left( \frac{-s}{\sqrt{2}} + \frac{1}{10} \right)
\left( \frac{1}{\sqrt{2}} - \frac{s}{10\,\sqrt{2}} \right)
\end{bmatrix} {\text d} s =
\begin{bmatrix}
\frac{\pi}{4} \\ \frac{1}{10} \\ \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} +
\begin{bmatrix} \frac{t}{10} - \frac{t^2}{2\,\sqrt{2}}
\\
-\frac{t}{\sqrt{2}} + \frac{t^2}{20\,\sqrt{2}}
\\
- \frac{t}{10\,\sqrt{2}} + \frac{t^2}{4} - \frac{t^2}{200\,\sqrt{2}} + \frac{t^3}{60}
\\
\frac{t}{10\,\sqrt{2}} - \frac{t^2}{4} - \frac{t^2}{200\,\sqrt{2}} + \frac{t^3}{60}
\end{bmatrix} .
\]
-(1/Sqrt[2] + t/10/Sqrt[2])*(1/10 - t/Sqrt[2])
-(1/(10 Sqrt[2])) + t/2 - t/(100 Sqrt[2]) + t^2/20
(1/Sqrt[2] - t/10/Sqrt[2])*(1/10 - t/Sqrt[2])
Expand[%]
Expand[%]
1/(10 Sqrt[2]) - t/2 - t/(100 Sqrt[2]) + t^2/20
We can expand this procedure indefinitely and obtain
\[
\theta (t) = \frac{\pi}{4} + \frac{t}{10} -\frac{\sqrt{2}}{4}\, t^2 - \frac{\sqrt{2}}{120}\, t^3 + \left( \frac{1}{48} + \frac{\sqrt{2}}{4800} \right) t^4 + \left( \frac{\sqrt{2}}{240 000} - \frac{1}{1200} \right) t^5 + \left( \frac{1111\,\sqrt{2}}{1600 000} - \frac{11}{144 000} \right) t^6 + \cdots .
\]
■
Example 4:
■
Example 5:
■
Example 6:
■
Example 7:
Consider the initial value problem for the homogeneous Duffing equation:
\[
\ddot{y} + 2\eta\,\dot{y} + y + \varepsilon \, y^3 = 0 , \qquad
y(0) = d, \quad \dot{y} (0) = v ,
\]
where d = 1 is the initial displacement and v =-1 is the initial
velocity. For simplicity, we choose the following numerical values of
parameters: η = 3 and ε = 1. First, we convert the given problem
to a system of first order differential equations
■
Example 8:
■
Example:
**DESCRIPTION OF PROBLEM GOES HERE**
This is a description for some MATLAB code. MATLAB is an extremely
useful tool for many different areas in engineering