一、簡介

常微分方程數值解法(numerical methods forordinary differential equations)計算數學的一個分支。是解常微分方程各類定解問題的數值方法，現有的解析方法只能用於求解一些特殊型別的定解問題，實用上許多很有價值的常微分方程的解不能用初等函式來表示，常常需要求其數值解。所謂數值解，是指在求解區間內一系列離散點處給出真解的近似值。

詳解

二、實現

# -*- coding: utf-8 -*-
"""
Created on Mon Dec 19 15:55:39 2016

    Euler 改進Euler Runge-Kutta

@author: Administrator
"""
 


from numpy import *
import matplotlib.pyplot as plt

#[0,1]
def f(t,u):
    return 1 / (t + 1)**2 * (t * u - u**2)

def euler(a,b,step):
    u = [] 
    t = []
    u.append(2)
    t.append(0)
    n  = int((b - a) / step)
    for i in range(0,n):
        t.append(a + (i + 1) * step)
        u.append(u[i] + step * f(t[i],u[i]))
    return
 
 t,u

def euler_ex(a,b,step):
    u = []
    t = []
    u.append(2)
    t.append(0)
    n  = int((b - a) / step)
    for i in range(0,n):
        t.append(a + (i + 1) * step)
        temp_u = u[i] + step * f(t[i],u[i])
        u.append(u[i] + step/2 * (f(t[i],u[i]) + f(t[i+1],temp_u))) 
    return t,u

#四階 Runge Kutta
 

def runge_kutta(a,b,step):
    u = []
    t = []
    u.append(2)
    t.append(0)
    n  = int((b - a) / step)
    for i in range(0,n):
        t.append(a + (i + 1) * step)
        k1 = f(t[i],u[i])
        k2 = f(t[i] + step / 2,u[i] + step * k1 / 2)
        k3 = f(t[i] + step / 2,u[i] + step * k2 / 2)
        k4 = f(t[i] + step,u[i] + step * k3)
        u.append(u[i] + step/6 * (k1 + 2*k2 + 2*k3 + k4))
    return t,u

if __name__ == '__main__':
    a = 0
    b = 1
    step = 0.01  #0.1 0.01 0.001
    euler_t,euler_u = euler(a,b,step)
    euler_ex_t,euler_ex_u = euler_ex(a,b,step)
    runge_kutta_t,runge_kutta_u = runge_kutta(a,b,step)

    fig = plt.figure(figsize=(8,4))
    ax = fig.add_subplot(111)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.plot(euler_t,euler_u,color='b',linestyle='-',label='Euler')
    ax.plot(euler_ex_t,euler_ex_u,color='r',linestyle='--',label='Euler EX')
    ax.plot(runge_kutta_t,runge_kutta_u,color='k',linestyle='-.',label='Runge-Kutta')
    ax.legend(loc='upper right')
    fig.show()
    fig.savefig('a.png')