常微分方程求解过程中ode45与C语言自编程计算结果比较

常微分方程求解过程中ode45与C语言自编程计算结果比较

---

前言

在求解常微分方程时，MATLAB中的ode求解器是必不可少的利器，几乎所有的常微分方程均可通过数值积分方法求解。其中最最常用的还要数ode45求解器，一般我们会先用ode45对我们要求解的方程试算一下，如果计算时间很长，再换用ode23等刚性求解器。众所周知，ode45的原理是四阶Runge-Kutta算法，有时由于我们需要得到更好的解，会不断的减小时间步长，所以有时ode求解器的求解速度不尽如人意。这篇文章对比了matlab中ode45求解器与自编程C语言Runge-Kutta算法的求解精度与求解速度。试验表明，C语言中Runge-Kutta算法计算的解与ode45解相差不大，但速度更快一些。

---

一、四阶Runge-Kutta算法原理

高阶常微分方程一般可以转化为一个一阶微分方程组。
例如：
$$\begin{gathered} y^{\prime \prime }=f(x,y,y^{\prime }) \\ y(x_0)=y_0,y^{\prime }(x_0)=y_0^{\prime } \end{gathered}$$
可以转化为
$$\begin{gathered} z=y^{\prime } \\ {z}^{\prime }=f(x,y,z) \\ 初始条件y(x_0)=y_0,z(x_0)=y_0^{\prime } \end{gathered}$$
Runge-Kutta格式为
$$\begin{gathered} y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4) \\ {z}_{n+1}=z_n+\frac{h}{6}(L_1+2L_2+2L_3+L_4) \end{gathered}$$
其中,
$$\begin{gathered} K_1=z_n \\ {L}_1=f(x_n,y_n,z_n) \\ {K}_2=z_n+\frac{h}{2}{L}_1 \\ {L}_2=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}{K}_1,z_n+\frac{h}{2}{L}_1) \\ {K}_3=z_n+\frac{h}{2}{L}_2 \\ {L}_3=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}{K}_2,z_n+\frac{h}{2}{L}_2) \\ {K}_4=z_n+h{L}_3 \\ {L}_4=f(x_{n+1},y_{n+1},z_{n+1}) \end{gathered}$$

二、C语言代码

代码如下（示例）：

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926
#define N 10000 //定义点数
#define _CRT_SECURE_NO_WARNINGS；
typedef float(*FUNTYPE)(float, float, float);
float f(float x, float y, float z)
{
	return float(-10*y+0.1);
}
float g(float x, float y, float z)
{
	return z;
}
float mul_runge_kutta_4(FUNTYPE f, FUNTYPE g, float sta, float end, float f0, float g0)
{
	float h = (end - sta) / N;
	float k1, k2, k3, k4, L1, L2, L3, L4;
	float x, y, z;
	FILE* filepath = fopen("E:\\filedata.txt", "w");
	int i;
	
	x = 0;
	y = f0;
	z = g0;
	for (i = 0; i <= N; i++)
	{
		x = i*h;  //x为一系列的点
		k1 = g(x, y, z);
		L1 = f(x, y, z);
		k2 = g(x + h / 2, y + h / 2 * k1, z + h / 2 * L1);
		L2 = f(x + h / 2, y + h / 2 * k1, z + h / 2 * L1);
		k3 = g(x + h / 2, y + h / 2 * k2, z + h / 2 * L2);
		L3 = f(x + h / 2, y + h / 2 * k2, z + h / 2 * L2);
		k4 = g(x + h, y + h*k3, z + h*L3);
		L4 = f(x + h, y + h*k3, z + h*L3);
		printf("x=%f,y=%f,z=%f,k1=%f\n", x, y, z,k1);
		fprintf(filepath, "%f %f %f\n",x, y, z);
		y = y + h*(k1 + 2 * k2 + 2 * k3 + k4) / 6;
		z = z + h*(L1 + 2 * L2 + 2 * L3 + L4) / 6;
		/***************消K后的格式*********************/
		/*L1=f(x,y,z);
		    L2=f(x+h/2,y+h/2*z,z+h/2*L1);
		    L3=f(x+h/2,y+h/2*z+h*h/4*L1,z+h/2*L2);
		    L4=f(x+h,y+h*z+h*h/2*L2,z+h/2*L3);
		y=y+h*z+h*h/6*(L1+L2+L3);
		   z=z+h/6*(L1+2*L2+2*L3+L4);
		/**********************************************/
	}
	fclose(filepath);
	return z;
}
int main(int argc, char *argv[])
{
	
	float a = 0, b = 30; //【区间下限，区间上限】
	int i; //循环计数
	printf("mul_runge_kutta=%f\n", mul_runge_kutta_4(f, g, a, b, 0.01, 0.01));
	getchar();
	return 0;
}

三、实例比较

实例1

$x^{\prime }=-w_n^{2}x+F$
这是一强迫振动方程，假设系统固有频率为$w_n^2=10$,扰动力$F=0.1$
从下图算计算结果可看到，上述C语言程序的计算结果与matlab中ode45计算结果相差不大，但是计算点是由我们等距给定，而ode45可以做到在函数变化剧烈的地方加密点。

实例2

$y_1^{{}^{\prime \prime }}-u(1-y_1^2)y_1^{{}^{\prime }}+y_1=0$其中,$u=1$.

两种方法计算结果依然相差不大。

---

总结

以上就是利用C语言自编程的计算效果与MATLAB中的ode45方法对比。两个实例表明了利用C语言自编的代码具有良好的精度，计算时间会更快，在某些简单的计算场合可以直接应用。