Skip to content

DA Usage Demos

Jump to bottom Edit New page
ZhJ Sun edited this page Sep 28, 2022 · 1 revision

Picard Iteration with DA

Picard iteration is a method to solve analytically nonlinear differential equations. The DA technique can be used to approximate the solutions of initial value problems in high orders by the Picard iteration scheme.

Picard iteration scheme for solving ordinary differential equations (ODE)

Consider the initial value problem $$\frac{dy}{dx}=f(x,y), \quad y(x_0)=y_0.$$

Suppose that the slope function $f$ satisfies the existence and uniqueness theorem (Picard--Lindelöf) in some rectangular domain $R=[x_0−a,x_0+a] \times [y_0−b,y_0+b]$:

  • The slope function is bounded: $M=\max_{(x,y) \in R} \left| f(x,y) \right|$.
  • The slope function is Lipschitz continuous: $|f(x,y)−f(x,z)|≤L|y−z|$, for all $(x,y)$ from $R$.

Then the successive approximation $$\phi_n(x)=y_0+\int^x_{x_0}f\big(t,\phi_{n−1}(t)\big)dt,\quad \phi_0=y_0,$$ converges to the unique solution of the initial value problem $y′=f(x,y), \ y(x_0)=y_0$, in the interval $x_0−h \le x \le x_0+h$, where $h=\min\big(a,\dfrac{b}{M}\big)$.

Example

/// This is an example implementation of Picard iteration method for solving ODE

#include <dace/dace.h>
#include <cmath>

using namespace std;
using namespace DACE;

// The first example ODE
DA fun_1( DA x, DA t) {
    // dx = x + t, x(0) = 0
    // the analytic solution is: x = e^t - t - 1
    return x + t;
}

// The second example ODE
AlgebraicVector<DA> fun_2( AlgebraicVector<DA> x, DA t ) {
    // dx = y, dy = t
    // the analytic solution is:
    // x = x0 + y0*t + 1/6*t^3
    // y = y0 + 1/2*t^2
    AlgebraicVector<DA> dx;
    dx.resize(2);
    dx.at(0) = x.at(1);
    dx.at(1) = t;

    return dx;
}

int main( void ) {
    DA::init( 12, 2 );

    // Example 1
    DA x0 = 0;
    DA t = DA(1);

    cout << "function: x = e^t - t - 1" << endl;
    DA xi, xj;
    xi = x0;
    for( int i = 0; i < 10; ++i ) {
        // iteration 10 times
        xj = x0 + fun_1(xi, t).integ(1);
        xi = xj;

        cout << "x_" << i+1 << endl << xj << endl;
    }

    // Example 2
    AlgebraicVector<DA> x_0(2);
    x_0.at(0) = 1;
    x_0.at(1) = 1;

    cout << "function: x = x0 + y0*t + 1/6*t^3, y = y0 + 1/2*t^2" << endl;
    AlgebraicVector<DA> x_i(2), x_j(2);
    x_i = x_0;
    for( int i = 0; i < 5; ++i ) {
        // iteration 5 times
        x_j = x_0 + fun_2(x_i, t).integ(1);
        x_i = x_j;

        cout << "iteration " << i+1 << endl;
        cout << "x" << endl << x_i.at(0) << endl;
        cout << "y" << endl << x_i.at(1) << endl;
    }

    return 0;
}

The running ouputs are:

function: x = e^t - t - 1
x_1
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
------------------------------------------------

x_2
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
------------------------------------------------

x_3
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
------------------------------------------------

x_4
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
------------------------------------------------

x_5
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
     5    1.3888888888888889e-03   6   6  0
------------------------------------------------

x_6
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
     5    1.3888888888888889e-03   6   6  0
     6    1.9841269841269841e-04   7   7  0
------------------------------------------------

x_7
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
     5    1.3888888888888889e-03   6   6  0
     6    1.9841269841269841e-04   7   7  0
     7    2.4801587301587302e-05   8   8  0
------------------------------------------------

x_8
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
     5    1.3888888888888889e-03   6   6  0
     6    1.9841269841269841e-04   7   7  0
     7    2.4801587301587302e-05   8   8  0
     8    2.7557319223985893e-06   9   9  0
------------------------------------------------

x_9
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
     5    1.3888888888888889e-03   6   6  0
     6    1.9841269841269841e-04   7   7  0
     7    2.4801587301587302e-05   8   8  0
     8    2.7557319223985893e-06   9   9  0
     9    2.7557319223985894e-07  10  10  0
------------------------------------------------

x_10
     I  COEFFICIENT              ORDER EXPONENTS
     1    5.0000000000000000e-01   2   2  0
     2    1.6666666666666666e-01   3   3  0
     3    4.1666666666666664e-02   4   4  0
     4    8.3333333333333332e-03   5   5  0
     5    1.3888888888888889e-03   6   6  0
     6    1.9841269841269841e-04   7   7  0
     7    2.4801587301587302e-05   8   8  0
     8    2.7557319223985893e-06   9   9  0
     9    2.7557319223985894e-07  10  10  0
    10    2.5052108385441720e-08  11  11  0
------------------------------------------------

function: x = x0 + y0*t + 1/6*t^3, y = y0 + 1/2*t^2
iteration 1
x
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   1  0
------------------------------------------------

y
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    5.0000000000000000e-01   2   2  0
------------------------------------------------

iteration 2
x
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   1  0
     3    1.6666666666666666e-01   3   3  0
------------------------------------------------

y
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    5.0000000000000000e-01   2   2  0
------------------------------------------------

iteration 3
x
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   1  0
     3    1.6666666666666666e-01   3   3  0
------------------------------------------------

y
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    5.0000000000000000e-01   2   2  0
------------------------------------------------

iteration 4
x
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   1  0
     3    1.6666666666666666e-01   3   3  0
------------------------------------------------

y
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    5.0000000000000000e-01   2   2  0
------------------------------------------------

iteration 5
x
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    1.0000000000000000e+00   1   1  0
     3    1.6666666666666666e-01   3   3  0
------------------------------------------------

y
     I  COEFFICIENT              ORDER EXPONENTS
     1    1.0000000000000000e+00   0   0  0
     2    5.0000000000000000e-01   2   2  0
------------------------------------------------
Add a custom footer
Add a custom sidebar
Clone this wiki locally