软件下载 | 资讯教程 | 最近更新 | 下载排行 | 一键转帖 | 发布投稿
您的位置:最火下载站 > 电脑教程 > 编程开发 > C/C++开发 > 数值计算:微分方程的数值解法的C++程序

数值计算:微分方程的数值解法的C++程序

微分方程的各种算法的C++实现,是学数值计算方法的实验,有欧拉法,二阶龙格库塔法和四阶龙格库塔法。 程序使用了我自己写的库(已上传)。

Linux下编译方法:g++ df.cpp -o df -lm

附图为把程序对y'=cos(x)在[0,10]区间内插值运行结果复制到scilab里画的图,实线是原始函数图像,点实线为欧拉法的解,*号为RK2的解,虚线为RK4结果。

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <mylib/InitData.h>

using namespace std;

//微分方程 df = f(x,y);
double df(double x, double y)
{
return cos(x); // df/dx = cos(x); 解为f(x) = sin(x) + C
}

//欧拉法解微分方程,分块数量为n, 解的区间为[a,b], 解向量为(x,y),方程初值为(x0,y0)
bool Eluer(vector<double> & x, vector<double> & y, const int & n,
const double & a, const double & b,
const double & x0, const double & y0)
{
double h = (b-a)/n;
x.clear();
y.clear();
x.resize(n+1);
y.resize(n+1);
x[0] = x0;
y[0] = y0;
for (int i = 1; i <= n; i++)
{
x[i] = a + h*i;
y[i] = x[i-1] + h*df(x[i-1], y[i-1]);
}
return true;
}

//二阶龙格库塔法解解微分方程,分块数量为n, 解的区间为[a,b], 解向量为(x,y),方程初值为(x0, y0)
//参考:http://blog.sina.com.cn/s/blog_698c6a6f0100lp4x.html
bool RK2(vector<double> & x, vector<double> & y, const int & n,
const double & a, const double & b,
const double & x0, const double & y0)
{
double h = (b-a)/n;
x.clear();
y.clear();
x.resize(n+1);
y.resize(n+1);
x[0] = x0;
y[0] = y0;
double K1 = 0.0, K2 = 0.0;
for (int i = 1; i <= n; i++)
{
x[i] = a + h*i;
K1 = df(x[i-1], y[i-1]);
K2 = df(x[i-1] + h, y[i-1] + h*K1);
y[i] = y[i-1] + h/2.0*(K1 + K2);
}
return true;
}

//四阶龙格库塔法解解微分方程,分块数量为n, 解的区间为[a,b], 解向量为(x,y),方程初值为(x0, y0)
//参考:http://blog.sina.com.cn/s/blog_698c6a6f0100lp4x.html 和维基百科
bool RK4(vector<double> & x, vector<double> & y, const int & n,
const double & a, const double & b,
const double & x0, const double & y0)
{
double h = (b-a)/n;
x.clear();
y.clear();
x.resize(n+1);
y.resize(n+1);
x[0] = x0;
y[0] = y0;
double K1 = 0.0, K2 = 0.0, K3 = 0.0, K4 = 0.0;
for (int i = 1; i <= n; i++)
{
x[i] = a + h*i;
K1 = df(x[i-1], y[i-1]);
K2 = df(x[i-1] + 0.5*h, y[i-1] + 0.5*h*K1);
K3 = df(x[i-1] + 0.5*h, y[i-1] + 0.5*h*K2);
K4 = df(x[i-1] + h, y[i-1] + h*K3);
y[i] = y[i-1] + h/6.0*(K1 + 2*K2 + 2*K3 + K4);
}
return true;
}

int main()
{
int n;
double a, b;
cout<<"输入求解区间:"<<endl;
cout<<"a = ";
cin>>a;
cout<<"b = ";
cin>>b;
cout<<"输入分块个数:";
cin>>n;
if (n <= 1 || a >= b)
{
cout<<"输入错误!"<<endl;
exit(1);
}
double x0, y0;
cout<<"输入方程初值:"<<endl;
cout<<"x0 = ";
cin>>x0;
cout<<"y0 = ";
cin>>y0;
vector<double> x(n+1), y(n+1); //解向量

//欧拉法
if (Eluer(x, y, n, a, b, x0, y0) == true)
{
cout<<"欧拉法求解成功:"<<endl;
cout<<"x = "<<endl;
ShowData(x, cout);
cout<<"y = "<<endl;
ShowData(y, cout);
}
else
{
cout<<"欧拉法求解失败"<<endl;
}

//二阶RK法
if (RK2(x, y, n, a, b, x0, y0) == true)
{
cout<<"二阶龙格库塔法求解成功:"<<endl;
cout<<"x = "<<endl;
ShowData(x, cout);
cout<<"y = "<<endl;
ShowData(y, cout);
}
else
{
cout<<"二阶龙格库塔法求解失败"<<endl;
}

//四阶RK法
if (RK4(x, y, n, a, b, x0, y0) == true)
{
cout<<"四阶龙格库塔法求解成功:"<<endl;
cout<<"x = "<<endl;
ShowData(x, cout);
cout<<"y = "<<endl;
ShowData(y, cout);
}
else
{
cout<<"四阶龙格库塔法求解失败"<<endl;
}
return 0;
}

    相关阅读
    栏目导航
    推荐软件