longe1.cpp

龙格库塔求解微分方程数值解

/*
/*编程者：刘艮平
/*完成日期：2004.5.29
/*e.g:y'=y-2x/y,x∈[0,0.6]
/*　　y(0)=1
/*使用经典四阶龙格-库塔算法进行高精度求解　
/* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6
/* K1=h*f(xi,yi)
/* K2=h*f(xi+h/2,yi+K1/2)
/* K3=h*f(xi+h/2,yi+K2/2)
/* K4=h*f(xi+h,yi+K3)
*/
#include <stdio.h>
#include <math.h>

float f(float den,float p0)  //要求解的微分方程的右部的函数 e.g: y-2x/y
{
	float rus;
 //   den=w/(W0+sl);
 //   rus=k*A*f/Wp*sqrt(RTp)-(k-1)*.....
	rus=p0-2*den/p0;
	return(rus);
}

//float fden()
//{

//}


void main()
{
	float x0; //范围上限
	int x1;   //范围下限
	float h;  //步长
	int n;    //计算出的点的个数
 	float k1,k2,k3,k4;
    float y[3]; //用于存放计算出的常微分方程数值解
	int i=0;
	int j;

//以下为函数的接口
	printf("intput x0:");
    scanf("%f",&x0);

	printf("input x1:");
	scanf("%f",&x1);

    printf("input y[0]:");
    scanf("%f",&y[0]);

	printf("input h:");
    scanf("%f",&h);

// 以下为核心程序
    n=((x1-x0)/h);
	n=3;
	
    for(j=0;j<n;j++)
	{   

      k1=h*f(x0,y[i]); //求K1
      k2=h*f((x0+h/2),(y[i]+k1/2)); //求K2
      k3=h*f((x0+h/2),(y[i]+k2/2)); //求K3
      k4=h*f((x0+h),(y[i]+k3)); //求K4
 
	  y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求yi+1
	  x0+=float(0.2);
      printf("y[%f]=%f\n",x0,y[i+1]); 
      ++i;
	 
	}
 
}