work.cpp

用vc实现数值分析中常微分初值问题的数值解法

// work.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "iostream.h"
#include "stdio.h"
#define f(x,y) (y*y)
#define g(x,y) (x/y)
//f(x,y),g(x,y)为所定义的函数表达式
#define N 20

int m;
double h;
double x[N],y[N],yb[N];
int main(int argc, char* argv[])
{	
	int i;
	int choice;
	double a,b,y0;
	void euler(double x[],double y[]);
	void rk(double x[],double y[]);
	cout<<"输入x的范围[a,b]，上下界a,b分别为："<<endl;
	cin>>a>>b;
	cout<<"y("<<a<<")=";
	cin>>y0;
	cout<<"输入将区间划分的快数(0<m<20)m=";
	cin>>m;
	while((m>N)||(m<0))
	{
		cout<<"所定义最大区间为20，请输入20以内整数m=";
		cin>>m;
	}
	h=(b-a)/m;
	x[0]=a;
	y[0]=y0;
	
	cout<<"请选择( 1 / 2 ) :1.改进的euler法,2.Runge-kutta法"<<endl;
	cin>>choice;
	if (choice==1)
	{
    	euler(x,y);
	   	for(i=0;i<=m;i++)
		{
			cout<<"x"<<i<<"="<<x[i]<<" , ";
		    cout<<"yb"<<i<<"="<<yb[i]<<" , ";
		    cout<<"y"<<i<<"="<<y[i]<<endl;
		}
	}
	else 
	{
		rk(x,y);
		for(i=0;i<=m;i++)
		{
			cout<<"x"<<i<<"="<<x[i]<<" , ";
		    cout<<"y"<<i<<"="<<y[i]<<endl;
		}
	}
	return 0;
}

	
//改进的euler法
void euler(double x[],double y[])
{  	   
	 int i;
	 for(i=1;i<=m;i++)
	 {
		x[i]=x[i-1]+h;
		yb[i]=y[i-1]+h*f(x[i-1],y[i-1]);
		y[i]=y[i-1]+h/2*(f(x[i-1],y[i-1])+f(x[i],yb[i]));
	 }
 }
//四阶Runge-Kutta法
void rk(double x[],double y[])
{
	int i;
	double k1,k2,k3,k4;
	for(i=1;i<=m;i++)
	{
		k1=h*g(x[i-1],y[i-1]);
		k2=h*g(x[i-1]+h/2,y[i-1]+k1/2);
		k3=h*g(x[i-1]+h/2,y[i-1]+k2/2);
		k4=h*g(x[i-1]+h,y[i-1]+k3);
		x[i]=x[i-1]+h;
		y[i]=y[i-1]+(k1+2*k2+2*k3+k4)/6;	
	}
}