# MATLAB和Python求解非线性常微分方程

## 龙格-库塔方法

M阶Runge-Kutta方法用于求解方程，

d y d x = f ( x , y ) , y ( a ) = y 0    ( 1 ) \frac{\boldsymbol{dy}}{\boldsymbol{dx}}=\boldsymbol{f}\left( \boldsymbol{x},\boldsymbol{y} \right) ,\boldsymbol{y}\left( \boldsymbol{a} \right) =\boldsymbol{y}_0\,\, \left( 1 \right) dxdy=f(x,y),y(a)=y0(1)

∑ j = 1 M a i j = c i \sum_{\boldsymbol{j}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{ij}}=\boldsymbol{c}_{\boldsymbol{i}}} j=1Maij=ci

∑ i = 1 M b i = 1 \sum_{\boldsymbol{i}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{i}}=1} i=1Mbi=1

y ( x i + 1 ) = y ( x i ) + ∫ x i x i + 1 f ( x , y ( x ) ) d x , i = 0 , ⋯   , N − 1 ( 2 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}+1} \right) =\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}} \right) +\int_{\boldsymbol{x}_{\boldsymbol{i}}}^{\boldsymbol{x}_{\boldsymbol{i}+1}}{\boldsymbol{f}\left( \boldsymbol{x},\boldsymbol{y}\left( \boldsymbol{x} \right) \right) \boldsymbol{dx}},\boldsymbol{i}=0,\cdots ,\boldsymbol{N}-1 \left( 2 \right) y(xi+1)=y(xi)+xixi+1f(x,y(x))dx,i=0,,N1(2)

f ( x i ( j ) , y ( x i ( j ) ) ) ≈ ∑ l = 1 M a j l f ( x i ( l ) , y ( x i ( l ) ) ) \boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) \approx \sum_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{jl}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)} \right) \right)} f(xi(j),y(xi(j)))l=1Majlf(xi(l),y(xi(l)))

y ( x i ( j ) ) = y ( x i ( 1 ) ) + ∫ x i ( 1 ) x i ( j ) f ( x , y ( x ) ) d x ≈ y ( x i ( 1 ) ) + h i ∑ l = 1 M a j l f ( x i ( l ) , y ( x i ( l ) ) )    ( 3 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) =\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)} \right) +\int_{\boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)}}^{\boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)}}{\boldsymbol{f}\left( \boldsymbol{x},\boldsymbol{y}\left( \boldsymbol{x} \right) \right) \boldsymbol{dx}}\\\approx \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)} \right) +\boldsymbol{h}_{\boldsymbol{i}}\sum_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{jl}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)} \right) \right)}\,\, \left( 3 \right) y(xi(j))=y(xi(1))+xi(1)xi(j)f(x,y(x))dxy(xi(1))+hil=1Majlf(xi(l),y(xi(l)))(3)

∑ l = 1 M b j f ( x i ( j ) , y ( x i ( j ) ) ) , ∑ i = 1 M b i = 1 \sum_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{j}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) ,\sum_{\boldsymbol{i}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{i}}}}=1 l=1Mbjf(xi(j),y(xi(j))),i=1Mbi=1

y ( x i ( M ) ) ≈ y ( x i ( 1 ) ) + h i ∑ j = 1 M b j f ( x i ( j ) , y ( x i ( j ) ) )    ( 4 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{M} \right)} \right) \approx \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)} \right) +\boldsymbol{h}_{\boldsymbol{i}}\sum_{\boldsymbol{j}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{j}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) \,\, \left( 4 \right)} y(xi(M))y(xi(1))+hij=1Mbjf(xi(j),y(xi(j)))(4)

ξ i = y ( x i ( l ) ) + h i ∑ j = 1 M b j f ( x i ( j ) , y ( t i ( j ) ) ) − y ( x i ( M ) ) \boldsymbol{\xi }_{\boldsymbol{i}}=\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)} \right) +\boldsymbol{h}_{\boldsymbol{i}}\sum\nolimits_{\boldsymbol{j}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{j}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{t}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) -\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{M} \right)} \right)} ξi=y(xi(l))+hij=1Mbjf(xi(j),y(ti(j)))y(xi(M))定义了局部截断误差M阶龙格库塔法。 如果M阶龙格库塔法的局部截断误差为 O ( h P + 1 ) \mathcal{O}\left( \boldsymbol{h}^{\boldsymbol{P}+1} \right) O(hP+1) ，其中 h = max ⁡ i = 0 , ⋯   , N − 1 { h i } \boldsymbol{h}=\underset{\boldsymbol{i}=0,\cdots ,\boldsymbol{N}-1}{\max}\left\{ \boldsymbol{h}_{\boldsymbol{i}} \right\} h=i=0,,N1max{ hi} ，则M级Runge-Kutta方法的应该为 P ∈ N \boldsymbol{P}\in \mathbb{N} PN阶 。

## Python解算器

原文作者：亚图跨际
原文地址: https://blog.csdn.net/jiyotin/article/details/113813205
本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。