四阶龙格库塔法求解微分方程
作者：PEZHANG
时间：2021.11.6
求解过程数学描述
四阶龙格库塔的求解过程可用如下数学公式描述：
k 1 = f ( t n , y n ) k_1=f\left( t_n,y_n \right) k1=f(tn,yn)
k 2 = f ( t n + h 2 , y n + h 2 k 1 ) k_2=f\left( t_n+\frac{h}{2},y_n+\frac{h}{2}k_1 \right) k2=f(tn+2h,yn+2hk1)
k 3 = f ( t n + h 2 , y n + h 2 k 2 ) k_3=f\left( t_n+\frac{h}{2},y_n+\frac{h}{2}k_2 \right) k3=f(tn+2h,yn+2hk2)
k 4 = f ( t n + h , y n + h k 3 ) k_4=f\left( t_n+h,y_n+hk_3 \right) k4=f(tn+h,yn+hk3)
y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) y_{n+1}=y_n+\frac{h}{6}\left( k_1+2k_2+2k_3+k_4 \right) yn+1=yn+6h(k1+2k2+2k3+k4)
例子
为验证程序的有效性，选取一个已知解析解的微分方程验证。
equ： y ′ = y y^{'}=y y′=y ，零初值状态，即 y ( 0 ) = 1 y(0)=1 y(0)=1。
可知解析解为： y = e x y=e^x y=ex
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('解析解','数值解');
结果分析
根据预期，算法应当逐渐收敛至稳态，但是实际的求解过程无法反应此过程。当求解区间变大时，出现算法出现轻微的发散现象，说明算法设计存在缺陷，原因尚未分析出来，后续理清思路后再补充。
参考文献
【1】https://www.jianshu.com/p/e4aa9a688959
【2】https://wenku.baidu.com/view/d69e8f1f77c66137ee06eff9aef8941ea76e4b2c.html
2022-10-26-四阶龙格库塔法计算程序【修正】
由于后续的工作需要使用数值计算方法，重写了四阶龙格库塔法，通过控制求解步长，可以有效的控制误差，上次遗留的发散问题仍未得到解决。
再次阅读之前写的程序，发现公布的算法是在反向验证龙格库塔算法的有效性，在解析解未知的前提下，算法无法进行正向求解。最新改进的算法已实现了正向求解。
%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('解析解','数值解');
title('h=0.01')
subplot(2,1,2);
error = exp(t)-y;
plot(t,error)
legend('误差');
通过对比上方的两张图片，可以发现，在求解步长设置为0.01时，求解误差已经非常小。
2022-12-29-四阶龙格库塔法C语言计算程序
sublime text的C语言环境配置
{"working_dir": "$file_path","cmd": "gcc -Wall \"$file_name\" -o \"$file_base_name\"","file_regex": "^(..[^:]*):([0-9]+):?([0-9]+)?:? (.*)$","selector": "source.c","variants":[{"name": "Run","shell_cmd": "gcc -Wall \"$file\" -o \"$file_base_name\" && start cmd /c \"${file_path}/${file_base_name} & pause\""}]
}例子1、已知： y ′ = y − 2 x y ， y ( 0 ) = 1 y^{'}=y-\frac{2x}{y}\text{，}y\left( 0 \right) =1 y′=y−y2x，y(0)=1
#include<stdio.h>
#include<math.h>#define NS 1double x[NS];
double y[NS];
void equ(double t, double *x, double *fx){//fx[0] = x[0]; //求解y'=yfx[0] = x[0] - 2*y[0] / x[0]; //求解y' = y-2x/y,y=1
}void RK4(double t, double *x, double hs){double fx[NS];double k1[NS], k2[NS], k3[NS], k4[NS], xk[NS];int i;equ(t, x, fx);for(i=0; i<NS; ++i){ k1[i] = fx[i] * hs;xk[i] = x[i] + k1[i]*0.5; //对应yy[i] = t+hs/2; //对应x}equ(t, xk, fx); // timer.t+hs/2., for(i=0; i<NS; ++i){ k2[i] = fx[i] * hs;xk[i] = x[i] + k2[i]*0.5;y[i] = t+hs/2;}equ(t, xk, fx); // timer.t+hs/2., for(i=0;i<NS;++i){ k3[i] = fx[i] * hs;xk[i] = x[i] + k3[i];y[i] = t+hs;}equ(t, xk, fx); // timer.t+hs, for(i=0; i<NS; ++i){ k4[i] = fx[i] * hs;x[i] = x[i] + (k1[i] + 2*(k2[i] + k3[i]) + k4[i])/6.0;}
}int main(){// double sumtime = 0.1;double hs = 0.2;double t = 0;int jj = 0;x[0] = 1;y[0] = 0;for (int i = 0; i < 5; ++i){RK4(t, x, hs);t = t + hs;jj = jj + 1;printf("第%d次循环：%f\n", jj, x[0]);}return 0;
}
代码执行结果
例子2、已知： y ′ = y ， y ( 0 ) = 1 y^{'}=y\text{，}y\left( 0 \right) =1 y′=y，y(0)=1
#include<stdio.h>
#include<math.h>
#include <cstdlib>
#define NS 1double x[NS];
//double y[NS];
void equ(double t, double *x, double *fx){fx[0] = x[0]; //求解y'=y,y(0)=1//fx[0] = x[0] - 2*y[0] / x[0];
}void RK4(double t, double *x, double hs){double fx[NS];double k1[NS], k2[NS], k3[NS], k4[NS], xk[NS];int i;equ(t, x, fx);for(i=0; i<NS; ++i){ k1[i] = fx[i] * hs;xk[i] = x[i] + k1[i]*0.5;}equ(t, xk, fx); // timer.t+hs/2., for(i=0; i<NS; ++i){ k2[i] = fx[i] * hs;xk[i] = x[i] + k2[i]*0.5;}equ(t, xk, fx); // timer.t+hs/2., for(i=0;i<NS;++i){ k3[i] = fx[i] * hs;xk[i] = x[i] + k3[i];}equ(t, xk, fx); // timer.t+hs, for(i=0; i<NS; ++i){ k4[i] = fx[i] * hs;x[i] = x[i] + (k1[i] + 2*(k2[i] + k3[i]) + k4[i])/6.0;}
}void write_data_to_file(FILE *fw);int main(){// double sumtime = 0.1;double hs = 0.01;double t = 0;int jj = 0;x[0] = 1;//y[0] =0;remove("RKdata.dat");FILE* fw;for (int i = 0; i < 300; ++i){t = t + 0.01;fw = fopen("RKdata.dat", "a");RK4(t, x, hs);fprintf(fw, "%f,%f,", t, x[0]);jj = jj + 1;printf("第%d次循环：%f\n", jj, x[0]);fclose(fw);}system("python ./datavis.py"); return 0;
}
python代码：datavis.py，python环境配置可以参考下方的参考文献【4】
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import matplotlib; matplotlib.use('TkAgg')c = pd.read_csv("RKdata.dat", sep = ",", header=None,encoding='utf-8')
c = c.iloc[:,0:len(c.columns)-1]
#print(c)
#df= c[[i%2==0 for i in range(len(c.index))]]
#tt = c.iloc[:,0:len(c.columns)-1]
t = c.iloc[:,c.columns%2 == 0]x = c.iloc[:,c.columns%2 == 1]t.to_csv('excel2txt.txt', sep='\t', index=False)
a = np.loadtxt('excel2txt.txt')x.to_csv('excel2txt1.txt', sep='\t', index=False)
b = np.loadtxt('excel2txt1.txt')#print(a.shape)
#print(b.shape)x = np.arange(0, 3, 0.01)
#print(x.shape)
y = np.exp(x)z = b[1,:]-y
plt.subplot(2, 1, 1)plt.plot(x,y,label='analytical solution')
plt.plot(a[1,:],b[1,:],label='numerical solution')
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(x,z,label='numerical error')
plt.legend()plt.savefig("test.jpg", dpi=300,format="jpg")
plt.show()
例子2代码执行结果
参考文献
【3】https://blog.csdn.net/Devid_/article/details/105855009
【4】https://zhuanlan.zhihu.com/p/64445558
【5】https://www.bilibili.com/video/BV1JQ4y1k77Z/?spm_id_from=333.999.0.0
