常微分方程求解过程中ode45与C语言自编程计算结果比较
前言
在求解常微分方程时,MATLAB中的ode求解器是必不可少的利器,几乎所有的常微分方程均可通过数值积分方法求解。其中最最常用的还要数ode45求解器,一般我们会先用ode45对我们要求解的方程试算一下,如果计算时间很长,再换用ode23等刚性求解器。众所周知,ode45的原理是四阶Runge-Kutta算法,有时由于我们需要得到更好的解,会不断的减小时间步长,所以有时ode求解器的求解速度不尽如人意。这篇文章对比了matlab中ode45求解器与自编程C语言Runge-Kutta算法的求解精度与求解速度。试验表明,C语言中Runge-Kutta算法计算的解与ode45解相差不大,但速度更快一些。
一、四阶Runge-Kutta算法原理
高阶常微分方程一般可以转化为一个一阶微分方程组。
例如:
y
′
′
=
f
(
x
,
y
,
y
′
)
y
(
x
0
)
=
y
0
,
y
′
(
x
0
)
=
y
0
′
y''=f(x,y,y')\\ y(x_0)=y_0,y'(x_0)=y'_0
y′′=f(x,y,y′)y(x0)=y0,y′(x0)=y0′
可以转化为
z
=
y
′
z
′
=
f
(
x
,
y
,
z
)
初
始
条
件
y
(
x
0
)
=
y
0
,
z
(
x
0
)
=
y
0
′
z=y' \\ z'=f(x,y,z)\\ 初始条件y(x_0)=y_0,z(x_0)=y'_0
z=y′z′=f(x,y,z)初始条件y(x0)=y0,z(x0)=y0′
Runge-Kutta格式为
y
n
+
1
=
y
n
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
z
n
+
1
=
z
n
+
h
6
(
L
1
+
2
L
2
+
2
L
3
+
L
4
)
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)\\
yn+1=yn+6h(K1+2K2+2K3+K4)zn+1=zn+6h(L1+2L2+2L3+L4)
其中,
K
1
=
z
n
L
1
=
f
(
x
n
,
y
n
,
z
n
)
K
2
=
z
n
+
h
2
L
1
L
2
=
f
(
x
n
+
1
2
,
y
n
+
h
2
K
1
,
z
n
+
h
2
L
1
)
K
3
=
z
n
+
h
2
L
2
L
3
=
f
(
x
n
+
1
2
,
y
n
+
h
2
K
2
,
z
n
+
h
2
L
2
)
K
4
=
z
n
+
h
L
3
L
4
=
f
(
x
n
+
1
,
y
n
+
1
,
z
n
+
1
)
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+hL_3\\L_4=f(x_{n+1},y_{n+1},z_{n+1})\\
K1=znL1=f(xn,yn,zn)K2=zn+2hL1L2=f(xn+21,yn+2hK1,zn+2hL1)K3=zn+2hL2L3=f(xn+21,yn+2hK2,zn+2hL2)K4=zn+hL3L4=f(xn+1,yn+1,zn+1)
二、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
′
=
−
w
n
2
x
+
F
x'=-w_n^2x+F
x′=−wn2x+F
这是一强迫振动方程,假设系统固有频率为
w
n
2
=
10
w_n^2=10
wn2=10,扰动力
F
=
0.1
F=0.1
F=0.1
从下图算计算结果可看到,上述C语言程序的计算结果与matlab中ode45计算结果相差不大,但是计算点是由我们等距给定,而ode45可以做到在函数变化剧烈的地方加密点。
实例2
y
1
′
′
−
u
(
1
−
y
1
2
)
y
1
′
+
y
1
=
0
y_1^{''}-u(1-y_1^2)y_1^{'}+y_1=0
y1′′−u(1−y12)y1′+y1=0其中,
u
=
1
u=1
u=1.
两种方法计算结果依然相差不大。
总结
以上就是利用C语言自编程的计算效果与MATLAB中的ode45方法对比。两个实例表明了利用C语言自编的代码具有良好的精度,计算时间会更快,在某些简单的计算场合可以直接应用。