Numerical calculation method chapter8 Numerical solutions of ordinary differential equations

• Numerical calculation method Chapter8. Numerical solutions of ordinary differential equations
  • 0. Problem description
  • 1. Euler The formula
    • 1. forward Euler The formula
    • 2. backward Euler The formula
    • 3. Trapezoidal formula
  • 2. Runge-Kutta Method
    • 1. Second order Runge-Kutta Method
    • 2. Higher order Runge-Kutta Method
      • 1. Third order Runge-Kutta Method
      • 2. Four order Runge-Kutta Method
    • 3. python Pseudo code implementation
  • 3. Linear multistep method
    • 1. The basic idea
    • 2. Adams The formula
  • 4. Numerical solution of ordinary differential equations
    • 1. Numerical solution of first order ordinary differential equations
    • 2. Numerical methods for higher order differential equations

0. Problem description

The problems examined in this chapter are described in the title , That is, the numerical solution of ordinary differential equations ：

$$\{\begin{array}{rl}\frac{dy}{dx}& =f(x,y)\\ y(x_0)& =y_0\end{array}$$

1. Euler The formula

1. forward Euler The formula

Euler Formula is a relatively direct way to solve the numerical solution of ordinary differential equations ：

$$\frac{dy}{dx}=\frac{\delta y}{\delta x}=f(x,y)$$

Thus there are ：

$$y+\delta y=f(x,y)\cdot \delta x$$

thus , We can approximate each by constantly solving $x$ Value under $y$ value .

$$y_{n+1}=y_n+h\cdot f(x_n,y_n)$$

give python The pseudo code is implemented as follows ：

def fwd_euler_fn(f, x0, y0, step=1e-3):
    def fn(x):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        x, y = x0, y0
        for _ in range(n):
            y += h * f(x, y)
            x += h
        return y
    return fn

2. backward Euler The formula

backward Euler Formula and forward Euler There is no essential difference between formulas , However, the fine-tuning formula is ：

$$y_{n+1}=y_n+h\cdot f(x_{n+1},y_{n+1})$$

But it's obviously calculating $y_{n+1}$ In fact, I can't know $f(x_{n+1},y_{n+1})$ Value , Therefore, we need to superimpose the idea of iteration ：

$$y_{n+1}^{(k+1)}=y_n^{(k+1)}+h\cdot f(x_{n+1},y_{n+1}^{(k)})$$

The above iterative formula is called Picard iteration .

Can prove that , When $h$ Enough hours ,Picard Iterative convergence .

give python The pseudocode is as follows ：

def bwd_euler_fn(f, x0, y0, step=1e-3):
    def fn(x, epsilon = 1e-6):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        xlist = [x0 + i*h for i in range(n)]
        ylist = [y0 for i in range(n)]
        for i in range(n-1):
            ylist[i+1] = ylist[i] + h * f(xlist[i], ylist[i])
        while True:
            zlist = deepcopy(ylist)
            for i in range(n-1):
                zlist[i+1] = zlist[i] + h * f(xlist[i+1], ylist[i+1])
            delta = abs(zlist[-1] - ylist[-1])
            ylist = zlist
            if delta < epsilon:
                break
        return zlist[-1]
    return fn

3. Trapezoidal formula

Trapezoidal formula is still based on differential difference quotient in essence , However, it is different from the direct use of differentiation before , The expression of integral is used more strictly here , namely ：

$$y_{n+1}=y_n+{\int }_{x_n}^{x_{n+1}}f(x,y)dx$$

then , Here the trapezoidal formula is used to approximate the integration process , Yes ：

$$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))$$

Here we will also refer to $y_{n+1}$ The solution of the problem , however , Different from backward Euler The pure iterative idea in the formula , Only one iteration is used here to approximate , namely ：

$$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_n+h\cdot f(x_n,y_n)))$$

Again , We give python The pseudo code is implemented as follows ：

def tapezoid_fn(f, x0, y0, step=1e-3):
    def fn(x):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        x, y = x0, y0
        for _ in range(n):
            y += h/2 * (f(x, y) + f(x+h, y+h*f(x, y)))
            x += h
        return y
    return fn

2. Runge-Kutta Method

1. Second order Runge-Kutta Method

Runge-Kutta Compared with the previous Euler Formula is a method with relatively higher accuracy .

Euler The formula comes from the definition of differential , and Runge-Kutta The method investigated Taylor an , Yes ：

$$\begin{array}{rl}{y}_{n+1}& =y_n+\sum _i\frac{h^i}{i!}{y}^{(i)}(x_n)\\ & =y_n+h\cdot y^{\prime }(x_n)+\frac{h^2}{2!}\cdot y^{\prime \prime }(x_n)+...\end{array}$$

We keep the second-order term, that is ：

$$\begin{array}{rl}{y}_{n+1}& =y_n+h\cdot y^{\prime }(x_n)+\frac{h^2}{2}{y}^{\prime \prime }(x_n)\\ & =y_n+h\cdot f(x_n,y_n)+\frac{h^2}{2}(\frac{\partial f}{\partial x}(x_n,y_n)+\frac{\partial f}{\partial y}(x_n,y_n)\cdot f(x_n,y_n))\end{array}$$

however , One problem here is that the two partial derivatives are not necessarily easy to solve , Therefore, the above expression cannot be called directly .

and Runge-Kutta The method is to use an approximate displacement formula to estimate it , So that there is no error between them in the second derivative range .

To be specific ：

$$\begin{array}{rl}\Delta & =y_{n+}-y_n\\ & =h\cdot y^{\prime }(x_n)+\frac{h^2}{2}{y}^{\prime \prime }(x_n)\\ & =h\cdot f(x_n,y_n)+\frac{h^2}{2}(\frac{\partial f}{\partial x}(x_n,y_n)+\frac{\partial f}{\partial y}(x_n,y_n)\cdot f(x_n,y_n))\\ & =c_1\cdot h\cdot f(x_n,y_n)+c_2\cdot h\cdot f(x_n+ah,y_n+bhf(x_n,y_n))\\ & =c_1\cdot h\cdot f(x_n,y_n)+c_2\cdot h\cdot [f(x_n,y_n)+\frac{\partial f}{\partial x}(x_n,y_n)\cdot ah+\frac{\partial f}{\partial y}(x_n,y_n)\cdot bhf(x_n,y_n)]\end{array}$$

Or say , Give parameters $c_1,c_2,a,b$ Make the following conditions meet ：

$$\{\begin{array}{rl}{y}_{n+1}& =y_n+h(c_1\cdot k_1+c_2\cdot k_2)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+ah,y_n+bh{k}_1)\end{array}$$

bring $y_{n+1}$ There is no error in the range of the second derivative .

or ：

$$\{\begin{array}{rl}(c_1+c_2)h\cdot f(x_n,y_n)& =h\\ c_{2}a{h}^2\cdot \frac{\partial f}{\partial x}(x_n,y_n)& =\frac{h^2}{2}\frac{\partial f}{\partial x}(x_n,y_n)\\ c_{2}b{h}^2\cdot \frac{\partial f}{\partial y}(x_n,y_n)& =\frac{h^2}{2}\frac{\partial f}{\partial y}(x_n,y_n)\end{array}$$

Simplify the above formula to get ：

$$\{\begin{array}{rl}{c}_1+c_2& =1\\ 2c_{2}a& =1\\ 2c_{2}b& =1\end{array}$$

As long as the above equations are satisfied , The corresponding parameter $a,b,c_1,c_2$ Both of them can make the above two formulas have no error within the second derivative range .

Give two groups of common second order Runge-Kutta The formula is as follows ：

1. $c_1=\frac{1}{2},c_2=\frac{1}{2},a=1,b=1$

   $$\{\begin{array}{rl}{k}_1& =f(x_n,y_n)\\ k_2& =f(x_n+h,y_n+h{k}_1)\\ y_{n+1}& =y_n+\frac{h}{2}(k_1+k_2)\end{array}$$

2. $c_1=0,c_2=1,a=\frac{1}{2},b=\frac{1}{2}$

   $$\{\begin{array}{rl}{k}_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1)\\ y_{n+1}& =y_n+h{k}_2\end{array}$$

2. Higher order Runge-Kutta Method

alike , We follow the above ideas , Give the general higher order Runge-Kutta The expression of the method is as follows ：

$$\{\begin{array}{rl}{y}_{n+1}& =y_n+h(c_1\cdot k_1+c_2\cdot k_2+...+c_m\cdot k_m)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+a_{1}h+y_n+b_{1,1}h\cdot k_1)\\ & ...\\ k_m& =f(x_n+a_{m-1}h+y_n+\sum _{i=1}^{m-1}{b}_{m-1,i}h\cdot k_i)\end{array}$$

Parameters $a,b,c$ You can make $y_n$ have $O(h^{m+1})$ Order accuracy .

or ,Taylor Expand until $h^m$ Both sides have the same expansion coefficient .

1. Third order Runge-Kutta Method

We give the third order Runge-Kutta Two sets of typical coefficients of the method are as follows ：

1. Coefficient combination （ One ）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+4k_2+k_3)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(x_n+h,y_n-h{k}_1+2h{k}_2)\end{array}$$

2. Coefficient combination （ Two ）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{4}(k_1+3k_3)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{k}_1)\\ k_3& =f(x_n+\frac{2}{3}h,y_n+\frac{2}{3}h{k}_2)\end{array}$$

3. Coefficient combination （ 3、 ... and ）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{9}(2k_1+3k_2+4k_3)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(x_n+\frac{3}{4}h,y_n+\frac{3}{4}h{k}_2)\end{array}$$

2. Four order Runge-Kutta Method

alike , We can give two groups of typical fourth order Runge-Kutta The formula is as follows ：

1. Coefficient combination （ One ）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_1)\\ k_3& =f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{k}_2)\\ k_4& =f(x_n+h,y_n+h{k}_3)\end{array}$$

2. Coefficient combination （ Two ）

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{8}(k_1+3k_2+3k_3+k_4)\\ k_1& =f(x_n,y_n)\\ k_2& =f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{k}_1)\\ k_3& =f(x_n+\frac{2}{3}h,y_n+\frac{1}{3}h{k}_1+h{k}_2)\\ k_4& =f(x_n+h,y_n+h{k}_1-h{k}_2+h{k}_3)\end{array}$$

3. python Pseudo code implementation

below , Let's give Runge-Kutta Methodical python The pseudo code is implemented as follows ：

def runge_kutta_fn(f, x0, y0, a, b, c, step=1e-3):
    def fn(x):
        h = step if x >= x0 else -step
        n = math.ceil((x - x0) / h)
        m = len(c)
        assert(len(a) == m-1 and all(len(b[i]) == i+1 for i in range(m-1)))
        k = [0 for _ in range(m)]
        x, y = x0, y0
        for _ in range(n):
            k[0] = f(x, y)
            for i in range(m-1):
                xi = x + a[i] * h
                yi = y
                for j in range(i+1):
                    yi += h * b[i][j] * k[j]
                k[i+1] = f(xi, yi)
            x += h
            y += sum(ci * ki * h for ci, ki in zip(c, k))
        return y
    return fn

3. Linear multistep method

1. The basic idea

The idea source of linear multistep method is still integral formula ：

$$y(x)=y(x_0)+{\int }_{x_0}^x{y}^{\prime }(t)dt$$

Before ,Euler Formula and Runge-Kutta The method is to estimate the value of the integral directly .

among ,Euler The formula is direct ${\int }_x^{x+h}{y}^{\prime }(t)dt=y^{\prime }(x)h$, and Runge-Kutta The method is to introduce a series of offsets to make $y(x)$ stay $n$ rank Taylor There is no error in expansion .

The approximate idea of linear multistep method is the idea of using the previous interpolation formula , Come on $y^{\prime }(x)$ To fit , Then use this fitting function to calculate the following integral value .

The description in the book is as follows ：

If you use the integral node $x_n,x_{n-1},...,x_{n-q}$ Construct interpolation polynomial approximation $y^{\prime }(x)$, In the interval $[x_{n-p},x_{n+1}]$ Calculate numerical integration on ${\int }_{x_{n-p}}^{x_{n+1}}{y}^{\prime }(x)dx$, Then it is called structural calculation $y_{n+1}$ The method is linear multistep .

Special ：

• If you use $x_n,x_{n-1},...,x_{n-q}$ Construct interpolation polynomials , Then the fitting function is called explicit Adams The formula ;
• If you use $x_{n+1},x_n,...,x_{n+1-q}$ Construct interpolation polynomials , Then the fitting function is called implicit Admas The formula ;

2. Adams The formula

here , Let's not deduce , Directly give some common Adams The formula is as follows ：

1. Second order explicit Adams The formula

   $$y_{n+1}=y_n+\frac{h}{2}(3f(x_n,y_n)-f(x_{n-1},y_{n-1}))$$

2. Second order implicit Adams The formula

   $$y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},y_{n+1}))$$

3. Third order display Adams The formula

   $$y_{n+1}=y_n+\frac{h}{12}(23f(x_n,y_n)-16f(x_{n-1},y_{n-1})+5f(x_{n-2},y_{n-2}))$$

4. Third order implicit Adams The formula

   $$y_{n+1}=y_n+\frac{h}{12}(5f(x_{n+1},y_{n+1}+8f(x_n,y_n)-f(x_{n-1},y_{n-1}))$$

5. Fourth order explicit Adams The formula

   $$y_{n+1}=y_n+\frac{h}{24}(55f(x_n,y_n)-59f(x_{n-1},y_{n-1})+37f(x_{n-2},y_{n-2})-9f(x_{n-3},y_{n-3}))$$

6. Fourth order implicit Adams The formula

   $$y_{n+1}=y_n+\frac{h}{24}(9f(x_{n+1},y_{n+1}+19f(x_n,y_n)-5f(x_{n-1},y_{n-1})+f(x_{n-2},y_{n-2}))$$

4. Numerical solution of ordinary differential equations

1. Numerical solution of first order ordinary differential equations

We give the expression of the initial value problem of the first order ordinary differential equation as follows ：

$$\{\begin{array}{rl}\frac{d{y}_1}{dx}& =f_1(x,y_1,y_2,...,y_m)\\ \frac{d{y}_2}{dx}& =f_2(x,y_1,y_2,...,y_m)\\ & ...\\ \frac{d{y}_m}{dx}& =f_m(x,y_1,y_2,...,y_m)\\ y_1(x_0)& ={\eta }_1\\ y_2(x_0)& ={\eta }_2\\ & ...\\ y_m(x_0)& ={\eta }_m\end{array}$$

In fact, the numerical solution of ordinary differential equations can be copied from the original model , There is no need to expand in detail .

Let's take the system of binary equations as an example , Give some common solutions ：

$$\{\begin{array}{rl}\frac{dy}{dx}& =f(x,y,z)\\ \frac{dz}{dx}& =g(x,y,z)\\ y(x_0)& =y_0\\ z(x_0)& =z_0\end{array}$$

1. Forward direction Euler The formula

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+hf(x_n,y_n,z_n)\\ z_{n+1}& =z_n+hg(x_n,y_n,z_n)\end{array}$$

2. Trapezoidal formula

   $$\{\begin{array}{rl}\widehat{y_{n+1}}& =y_n+hf(x_n,y_n,z_n)\\ \widehat{z_{n+1}}& =z_n+hg(x_n,y_n,z_n)\\ y_{n+1}& =y_n+\frac{h}{2}(f(x_n,y_n,z_n)+f(x_{n+1},\widehat{y_{n+1}},\widehat{z_{n+1}}))\\ z_{n+1}& =z_n+\frac{h}{2}(g(x_n,y_n,z_n)+g(x_{n+1},\widehat{y_{n+1}},\widehat{z_{n+1}}))\end{array}$$

3. Four order Runge-Kutta Method

   $$\{\begin{array}{rl}{y}_{n+1}& =y_n+\frac{h}{6}(k_{11}+2k_{12}+2k_{13}+k_{14})\\ z_{n+1}& =z_n+\frac{h}{6}(k_{21}+2k_{22}+2k_{23}+k_{24})\\ k_{11}& =f(x_n,y_n,z_n)\\ k_{21}& =g(x_n,y_n,z_n)\\ k_{12}& =f(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{11},z_n+\frac{h}{2}{k}_{21})\\ k_{22}& =g(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{11},z_n+\frac{h}{2}{k}_{21})\\ k_{13}& =f(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{12},z_n+\frac{h}{2}{k}_{22})\\ k_{23}& =g(x_n+\frac{h}{2},y_n+\frac{h}{2}{k}_{12},z_n+\frac{h}{2}{k}_{22})\\ k_{14}& =f(x_n+h,y_n+h{k}_{13},z_n+h{k}_{23})\\ k_{24}& =g(x_n+h,y_n+h{k}_{13},z_n+h{k}_{23})\end{array}$$

2. Numerical methods for higher order differential equations

here , Let's look at the numerical solution of higher-order differential equations in one variable .

$$\{\begin{array}{rl}\frac{d^{m}y(x)}{d{x}^m}& =f(x,y,y^{\prime },...,y^{(m-1)})\\ y(x_0)& ={\eta }_0\\ y^{\prime }(x_0)& ={\eta }_1\\ & ...\\ y^{(m-1)}(x_0)& ={\eta }_{m-1}\end{array}$$

In fact, this kind of problem can be used as an application example of the above first-order ordinary differential equations , We only need to make the following transformations to completely transform the problem into a first-order ordinary differential equation system , Then we can use the previous numerical solution of first-order ordinary differential equations to solve .

$$\{\begin{array}{rl}\frac{d{y}_0}{dx}& =y_1(x)\\ \frac{d{y}_1}{dx}& =y_2(x)\\ & ...\\ \frac{d{y}_{m-2}}{dx}& =y_{m-1}(x)\\ \frac{d{y}_{m-1}}{dx}& =f(x,y_0,y_1,...,y_{m-1})\\ y_0(x_0)& ={\eta }_0\\ y_1(x_0)& ={\eta }_1\\ & ...\\ y_{m-1}(x_0)& ={\eta }_{m-1}\end{array}$$