四阶龙格库塔法求解微分方程【MATLAB||C】

四阶龙格库塔法求解微分方程

---

作者：PEZHANG
时间：2021.11.6

---

求解过程数学描述

四阶龙格库塔的求解过程可用如下数学公式描述：
$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)$

例子

为验证程序的有效性，选取一个已知解析解的微分方程验证。

​ equ：$y^{{}^{\prime }}=y$ ，零初值状态，即$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:N
    k1 = 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)=1

clc;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-1
    k1 = 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^{{}^{\prime }}=y-\frac{2x}{y}\text{，}y\left(0\right)=1$

#include<stdio.h>
#include<math.h>

#define NS 1

double x[NS];
double y[NS];
void equ(double t, double *x, double *fx){
	//fx[0] = x[0];   //求解y'=y
	fx[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; //对应y
        y[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^{{}^{\prime }}=y\text{，}y\left(0\right)=1$

#include<stdio.h>
#include<math.h>
#include <cstdlib>
#define NS 1


double 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