常微分方程数值解法

一、简介

常微分方程数值解法(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')