本次求解微分方程,用4阶龙格库塔方法,该方法的代数精度高,比欧拉法,三阶龙格库塔都高,本次的实现很粗糙,由于本人没有熟练掌握函数指针的方式,对于不同微分方程,并不能做到一个普适的输入版本,所以本算法只能对本次的微分方程适用,等掌握了函数指针的用法,再来重写,这次的算法很简陋,也没有做一些越界的异常处理、判断。
代码如下:
//龙哥库塔方法求解微分方程
/*
*4阶龙格库塔方法,运用了经典公式
*微分方程如下: y'=y-2x/y;
y(0)=1;
*函数的输入:数组x是自变量区间数组,数组y是自变量区间y,h为步长的大小,begin为区间的起始,end为区间的末端
*进行测试的时候,需要自己估算区间的大小,暂时还没有实现自动的功能~.~
*/
#include<iostream>
using namespace std;
void Runge_Kutta(double *x, double* y,double h,double begin,double end)
{
double K1, K2, K3, K4;
for (int i = 0; i <= (end - begin) / h; i++)
{
x[i] = begin + h*i;
K1 = y[i] - 2 * x[i] / y[i];
cout << K1 << " ";
K2 = y[i] + h / 2 * K1 - 2 * (x[i] + h / 2) / (y[i] + h / 2 * K1);
cout << K2 << " ";
K3 = y[i] + h / 2 * K2 - 2 * (x[i] + h / 2) / (y[i] + h / 2 * K2);
cout << K3 << " ";
K4 = y[i] + h*K3 - 2 * (x[i] + h) / (y[i] + h*K3);
cout << K4 << " ";
y[i + 1] = y[i] + h / 6 * (K1 + 2 * K2 + 2 * K3 + K4);
cout << y[i+1]<<endl << endl;
}
}
int main()
{
double a[6];
a[0] = 0.0;
double y[7];
y[0] = 1;
Runge_Kutta(a, y, 0.2, 0.0, 1.0);
for (int i = 0; i <7; i++)
cout << y[i] << endl;
}