欧拉法求解常微分方程(c++)【转载】

 摘自《c++和面向对象数值计算》,代码简洁明快,采用类进行封装实现代码,增强代码的重用性,通过继承可实现代码的重用,采用函数指针,通用性增强,在函数改变时只需要单独改变函数部分的代码,无需对所有代码进行重写,对其中代码稍加改动,源代码有缺陷,没有释放内存,会造成内存泄露,在构造函数当中,将源代码的在函数体中赋值,改为列表赋值,可以稍微提高代码的执行效率。 

#include<iostream>

#include <cmath>

#include <algorithm>

using namespace std;

   

class ode

{

      double tini;      //初始时间

      double ison;      //初始解

      double tend;      //结束时间

      double(*sfn)(double t, double x);      //源函数,此处采用函数指针的方式,更具通用性

public:

      ode(double t0, double x0, double T, double(*f)(double, double)) :      //类的构造函数

            tini(t0), ison(x0), tend(T), sfn(f){}

      double* euler(int n) const;            //n子区间欧拉方法

};

   

double f(double t, double x)      //源函数

{

      return x*(1 - exp(t)) / (1 + exp(t));

}

   

double exact(double t)      //真实解

{

      return 12 * exp(t) / pow(1 + exp(t), 2);

}

   

int main()

{

      ode exmp(0, 3, 2, f);      //x(0)=3,T=2

      double* soln = exmp.euler(100);            //100个子区间

      double norm = 0;

      double h = 2.0 / 100;      //计算步长

      for (int k = 1; k <= 100; k++)

            norm = max(norm, fabs(exact(k*h) - soln[k]));

      cout << "max norm of error by euler's method = "

            << norm << endl;

      delete[] soln;

      return 0;

}

   

double* ode::euler(int n) const            //显式欧拉方法

{

      double* x = new double[n + 1];      //近似解

      double h = (tend - tini) / n;      //步长

      x[0] = ison;            //初始解

      for (int k = 0; k < n; k++)

            x[k + 1] = x[k] + h*sfn(tini + k*h, x[k]);

      return x;

}
 

posted @ 2015-12-04 02:58  硫酸亚铜  阅读(1129)  评论(0编辑  收藏  举报