常微分方程的数值解-欧拉、四阶龙格-库塔法等C语言

迈不过友情╰ 2022-11-21 06:50 235阅读 0赞

\*题目描述：\*求常微分方程y’=y-2x/y && y(0)=1的数值解。解析解为：y=sqrt(1+2x)  
***编译环境vc++6.0***：  
代码如下：

/* 1、y'=y-2x/y && y(0)=1; 2、解析解为y=sqrt(1+2x) */
    
    #include <stdio.h>
    #include <math.h>
    #define N 5 //离散点个数,分别取5,8,10，对应步长0.2,0.125,0.1
    #define low 0 //定义域x区间下限
    #define high 1.0 //定义域x区间上限
    #define h (high-low)/N //步长
    
    /*************************4、解析解实现**********************/
    double AnalyticalSolution(double x) { 
    	return sqrt(1 + 2 * x);
    }
    
    /*************************1、显示欧拉实现**********************/
    double EulerMethod(double x)
    { 
    	//若x==0
    	if (fabs(x - low) < 1e-6)  //浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
    		return 1;				//f(low)== f(0)==1
    	else return EulerMethod(x - h) + h * (EulerMethod(x - h) - (2 * (x - h) / EulerMethod(x - h)));
    }
    
    /*************************2、隐式欧拉实现**********************/
    //隐式欧拉实现由下面两个函数组成
    double ImplicitEulerMethod_1(double x)			/*隐示欧拉第一次迭代, 用显示欧拉计算f(x)的初值，带入隐式欧拉方程 f(x)=f(x-1)+h*(f(x)-2*x/f(x))*/
    { 
    	if (fabs(x - low) < 1e-6)	//浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
    		return 1;
    	else
    		return ImplicitEulerMethod_1(x - h) + h * (EulerMethod(x) - 2 * x / EulerMethod(x));   //迭代
    }
    
    double ImplicitEulerMethod_2(double x)	/*隐示欧拉第2--n次迭代， 用第一次迭代的结果f(x)作为初值带入隐式欧拉方程进行第二次迭代， 第二次迭代的结果f(x)作为初值带入隐式欧拉方程进行第三次迭代..... 直到本次迭代的结果y1与上次的迭代结果y0，满足|y1-y0|<1e-6,则认为收敛 迭代结束 */
    { 
    	if (fabs(x - low) < 1e-6)  //浮点数之间不能直接判断是否相等，则可以用 fabs(f1-f2)<=给定差值精度 来判断f1和f2是否相等
    		return 1;
    	else { 
    		double y0 = ImplicitEulerMethod_1(x), y1 = ImplicitEulerMethod_1(x);		//将第一次迭代的结果作为初值带入第二次迭代
    		do { 
    			y0 = y1;
    			y1 = ImplicitEulerMethod_2(x - h) + h * (y0 - 2 * x / y0);
    		} while (fabs(y1 - y0) < 1e-6);		//本次迭代结果减去上一次迭代结果fabs(y1-y0)<=1e-6 认为是收敛
    		return y1;
    	}
    }
    /*************************3、两步欧拉实现**********************/
    double DoubleStepEuler(double x)			//两步欧拉
    { 
    	if (fabs(x - low) < 1e-6)	//若x==0 第一个初值给定
    		return 1;
    	else if (fabs(x - low - h) < 1e-6)   //若x==low+h 第二个初值用显示欧拉求出
    		return EulerMethod(x);
    	else return DoubleStepEuler(x - 2 * h) + 2 * h * (DoubleStepEuler(x - h) - 2 * (x - h) / DoubleStepEuler(x - h));
    }
    
    
    
    /*************************四阶龙格-库塔法**********************/
    double C_R_K_M(double x)
    { 
    	if (fabs(x - low) < 1e-6)	//若x==0 第一个初值给定
    		return 1;
    	else { 
    		double k1 = C_R_K_M(x - h) - 2 * (x - h) / C_R_K_M(x - h);
    		double k2 = C_R_K_M(x - h) + h / 2.0 * k1 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h / 2.0 * k1);
    		double k3 = C_R_K_M(x - h) + h / 2.0 * k2 - (2 * (x - h) + h) / (C_R_K_M(x - h) + h / 2.0 * k2);
    		double k4 = C_R_K_M(x - h) + h * k3  - (2 * (x - h) + h) / (C_R_K_M(x - h) + h  * k3);
    		return C_R_K_M(x - h) + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4);
    	}
    }
    
    
    
    /*************************4、格式化输出函数**********************/
    void prit(double (*p)(double x))      
    { 
    	int count=0;
    	double i;
    	for (i = low; high - i >0 ; i += h) { 
    		printf("y(%0.3lf)=%0.5lf\t", i, p(i));
    		count++;
    		if (count % 4 == 0)
    			printf("\n");
    
    	}
    	printf("y(%0.3lf)=%0.5lf\t", i, p(i));		//high==i时
    	printf("\n\n\n");
    }
    
    
    int main()
    { 
    	printf("解析解格式化输出:\n");
    	prit(AnalyticalSolution);
    	
    	printf("显示欧拉法格式化输出:\n");
    	prit(EulerMethod);
    
    	printf("隐示欧拉法格式化输出:\n");
    	prit(ImplicitEulerMethod_2);
    	
    	printf("两步欧拉法格式化输出:\n");
    	prit(DoubleStepEuler);
    	
    	printf("四阶龙格-库塔法格式化输出:\n");
    	prit(C_R_K_M);
    	return 0;
    }

输出结果：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zOTc5MTY2NQ_size_16_color_FFFFFF_t_70_pic_center]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl8zOTc5MTY2NQ_size_16_color_FFFFFF_t_70_pic_center]: /images/20221120/beebe639d15e4e19ad9d4d2a962fc08d.png