数值微分.cpp

来自「这是一些c++例程」· C++ 代码 · 共 65 行

// 数值微分.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#include "iostream.h"

class CDiff
{
public:
	double h, y;						//成员变量，分别表示运算步长和计算结果
	double k1, k2, k3, k4;				//成员变量，表示公式运算中间结果。
	double f(double, double);				//成员函数，描述微分方程的表达式
	double Runge();			 		//成员函数，进行龙格-库塔公式的计算，并输出最终结果。
};

double CDiff::f(double x, double y)			//成员函数，描述微分方程的表达式
{
	return -20*y+20*sin(x)+cos(x);
}

double CDiff::Runge()					//成员函数，进行龙格-库塔公式的计算，并输出最终结果。
{
		double y=1, h=0.01;					//y(0)初值为1，h第一段步长为0.01
		double k1,k2,k3,k4;
		for (double t=0;t<0.2;t+=h)
		{
			k1=f(t,y);
			k2=f(t+h/2,y+h*k1/2);
			k3=f(t+h/2,y+h*k2/2);
			k4=f(t+h,y+h*k3);
			y+=h/6*(k1+2*k2+2*k3+k4);
		};
		h=0.075;						//h第二段步长为0.075
		for (t=0.2;t<=1-h;t+=h)
		{
			k1=f(t,y);
			k2=f(t+h/2,y+h*k1/2);
			k3=f(t+h/2,y+h*k2/2);
			k4=f(t+h,y+h*k3);
			y+=h/6*(k1+2*k2+2*k3+k4);
		};
		h=1-t;
		k1=f(t,y);
		k2=f(t+h/2,y+h*k1/2);
		k3=f(t+h/2,y+h*k2/2);
		k4=f(t+h,y+h*k3);
		y+=h/6*(k1+2*k2+2*k3+k4);
		//返回y，为最终积分数值
	return y;
}

void main()
{
	CDiff test;						     //生成CDiff类的对象test
		double res=test.Runge();				//计算积分
		printf("y(1)= %.16f\n", res);			//输出结果
		double e=exp(-20.0)+sin(1.0)-res;
		printf("e= %.16f\n", exp(-20.0)+sin(1.0)-res);	//输出误差
		cout<<100*e/res<<endl;
}