数值计算方法实验四

1.用欧拉法、改进欧拉法和四阶龙格-库塔法求解一阶常微分方程初值问题

令步长h=0.1。

要求：

(1) 步长h的值从键盘终端输入；

(2) 将三种方法计算的每一个节点的函数值打印输出，并和精确值进行比较（精确解对应的解析表达式为），根据计算结果分析这三种方法求解结果的精度差别。

#include<iostream>
using namespace std;
int main(){
	double h,Y=1;
	int n;
	cout<<"输入步长h: "<<endl;
	scanf("%lf",&h);
	cout<<"欧拉法："<<endl;
	cout<<"  xn  "<<"    yn  "<<endl; 
	for(n=1;n<=10;n++){
	   Y = Y + h*(Y - 2*((n-1)*h)/Y);
	   printf("  %.1lf   %.6lf\n",n*0.1,Y);	
	}
	cout<<"改进欧拉法："<<endl; 
	double x0=0,y0=1,x1,y1,yp,yc;
	n=0;
	cout<<"  xn  "<<"    yn  "<<endl;
	do{
		x1 = x0 + h;
		yp = y0 + h*(y0 - 2*x0/y0);
		yc = y0 + h*(yp - 2*x1/yp);
		y1 = (yp+yc)/2;	
        printf("  %.1lf   %.6lf\n",x1,y1);
		n++;
		x0 = x1;
		y0 = y1;
	}while(n<10);
	x0=0,y0=1;
	n=0;
	double k1,k2,k3,k4;
	cout<<"四阶龙格-库塔法："<<endl; 
	cout<<"  xn  "<<"    yn  "<<endl;
	do{
		k1 = y0 - 2*x0/y0;
		k2 = y0 + h/2*k1 - (2*x0+h)/(y0+h/2*k1);
		k3 = y0 + h/2*k2 - (2*x0+h)/(y0+h/2*k2);
		k4 = y0 + h*k3   - 2*(x0+h)/(y0+h*k3);
		y1 = y0 + (h/6)*(k1+2*k2+2*k3+k4); 
		x1 = x0 + h;
		printf("  %.1lf   %.6lf\n",x1,y1);
		n++;
		x0 = x1;
		y0 = y1;
	}while(n<10);
	return 0;
}

2.用四阶龙格-库塔法求解一阶常微分方程组初值问题

x∈[0, 2]，h=0.1。

要求：

(1) 步长h的值从键盘终端输入；

(2) 将每一步的计算结果打印输出。

程序源代码：
#include<iostream>
using namespace std;
double f1(double x,double y1,double y2){
	return (1/(y2 - x));
}
double f2(double x,double y1,double y2){
	return (-1/y1 + 1);
} 
int main(){
	double h,x0=0,x1,y0=1,y1,y20=1,y21,k1,k2,k3,k4,l1,l2,l3,l4;
    int n=0;
    cout<<"请输入步长h:"<<endl;
	cin>>h;
	cout<<endl<<"每一步的计算结果为:"<<endl;
	while(n!=(2/h)){
		x1 = x0 + h;
		k1 = f1(x0,y0,y20);
		l1 = f2(x0,y0,y20);
		k2 = f1(x0+h/2,y0+h*k1/2,y20+h*l1/2);
		l2 = f2(x0+h/2,y0+h*k1/2,y20+h*l1/2);
		k3 = f1(x0+h/2,y0+h*k2/2,y20+h*l2/2);
		l3 = f2(x0+h/2,y0+h*k2/2,y20+h*l2/2);
		k4 = f1(x1,y0+h*k3,y20+h*l3);
		l4 = f2(x1,y0+h*k3,y20+h*l3);
		y1 = y0 + h*(k1+2*k2+2*k3+k4)/6;
		y21 = y20 + h*(l1+2*l2+2*l3+l4)/6;
		printf("x%d = %.6lf   y1%d = %.6lf   y2%d = %.6lf\n",n+1,x1,n+1,y1,n+1,y21);
		n++;
		x0 = x1;
		y0 = y1;
		y20 = y21;
	}
	return 0;
}