四階龍格庫塔法求解微分方程

求解過程數(shù)學描述

四階龍格庫塔的求解過程可用如下數(shù)學公式描述：
$k_1=f\left(t_n,y_n\right)$

$k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_1\right)$

$k_3=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}{k}_2\right)$

$k_4=f\left(t_n+h,y_n+h{k}_3\right)$

$y_{n+1}=y_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)$

例子

為驗證程序的有效性，選取一個已知解析解的微分方程驗證。

? $y^{{}^{\prime }}=y$ ，零初值狀態(tài)，即$y(0)=1$，$y=e^x$

code

編寫的MATLAB程序如下：

% function RK4() clc;clear; Ts = 0.01; h = Ts; time = 1.5; N = time/Ts; t = linspace(Ts,time,N);y = zeros(1,N+1); for m=2:Nk1 = exp(m*Ts);k2 = exp(m*Ts+h/2*k1);k3 = exp(m*Ts+h/2*k2);k4 = exp(m*Ts+h*k3);y(1,m+1) = y(1,m) +(k1+2*k2+2*k3+k4)*h/6; end figure plot(t,exp(t)) hold on y = y(1,2:N+1); y = y+1; plot(t,y) legend('解析解','數(shù)值解');

結果分析

根據(jù)預期，算法應當逐漸收斂至穩(wěn)態(tài)，但是實際的求解過程無法反應此過程。當求解區(qū)間變大時，出現(xiàn)算法出現(xiàn)輕微的發(fā)散現(xiàn)象，說明算法設計存在缺陷，原因尚未分析出來，后續(xù)理清思路后再補充。

參考文獻

https://www.jianshu.com/p/e4aa9a688959

https://wenku.baidu.com/view/d69e8f1f77c66137ee06eff9aef8941ea76e4b2c.html

---

2022-10-26-四階龍格庫塔法計算程序修正

由于后續(xù)的工作需要使用數(shù)值計算方法，重寫了四階龍格庫塔法，通過控制求解步長，可以有效的控制誤差，上次遺留的發(fā)散問題仍未得到解決。

再次閱讀之前寫的程序，發(fā)現(xiàn)公布的算法是在反向驗證龍格庫塔算法的有效性，在解析解未知的前提下，算法無法進行正向求解。最新改進的算法已實現(xiàn)了正向求解。

%eqution %y'=y y(0)=1clc;clear;% set the solution range and solution step h = 0.01; time = 5; N = time/h; t = linspace(h,time,N); y = zeros(1,N); y(1,1) = 1;% RK4 for m=1:N-1k1 = h*y(1,m);k2 = h*(y(1,m)+k1/2);k3 = h*(y(1,m)+k2/2);k4 = h*(y(1,m)+k3);y(1,m+1) = y(1,m) +(k1+2*k2+2*k3+k4)/6; end% data visualization figure subplot(2,1,1); plot(t,exp(t)) hold on plot(t,y) legend('解析解','數(shù)值解'); title('h=0.01') subplot(2,1,2); error = exp(t)-y; plot(t,error) legend('誤差');

通過對比上方的兩張圖片，可以發(fā)現(xiàn)，在求解步長設置為0.01時，求解誤差已經非常小。

總結

以上是生活随笔為你收集整理的四阶龙格库塔法求解微分方程【MATLAB】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。