摘自《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;
}