常微分方程数值解法课程设计.doc
课程设计报告课程设计题目： 常微分方程数值解法学生姓名: 专 业：班 级： 指导教师： 时 间：题目：常微分方程数值解法用欧拉方法、改进欧拉方法、3阶龙格库塔法以及4阶龙格库塔法求解常微分方程初值问题：一、摘要在matlab环境下熟悉的运用计算机编程语言并结合龙格库塔法的理论基础对常微分方程初值问题进行求解，在运行完程序后以及对运行结果做出各方面的分析和比较。二、设计目的用熟悉的计算机语言编程上机完成用欧拉方法、改进欧拉方法、3阶龙格库塔法以及4阶龙格库塔法求解常微分方程初值问题。二、理论基础1.欧拉公式【3】在点 将作Taylor展开，得，那么当h充分小时，略去误差项，用近似替代、近似替代，并注意到，便得上述方法称为Euler方法。2.改进Euler方法【3】在应用梯形方法的迭代公式进行运算时，每迭代一次都要重新计算函数的值，且还要判断何时可以终止或转下一步计算。为了控制计算量和简化算法，通常只迭代一次就转入下一步计算。具体说，我们先用Euler公式求得一个初步的近似值，称之为预测值，然后用梯形方法的迭代公式作一次迭代得，即将校正一次，这样建立的预测-校正方法称之为改进的Euler方法：预测：，校正：.3.三阶龙格-库塔方法【3】类似前面改进的Euler方法公式的推导方法，将在处作Taylor展开，然后再将在处作Taylor展开，只要将两个展开式前四项相同便有。于是得到三阶龙格-库塔公式为：4.四阶龙格-库塔方法【3】 类似前面三阶龙格-库塔的推导方法，如果每步计算四次函数f（x，y）的值，完全类似的，可以导出局部截断误差为的四阶龙格-库塔公式，其公式为：三、程序代码及运算结果用欧拉方法、改进欧拉方法、3阶龙格库塔法以及4阶龙格库塔法求解常微分方程初值问题：1. 求精确解首先可以求得其精确解为：y=exp(-cos(x)>> x=0:0.05:1;>> y=exp(-cos(x);>> plot(x,y,'b*-');>> Data=x',y' Data = 0 0.3679 0.0500 0.3683 0.1000 0.3697 0.1500 0.3720 0.2000 0.3753 0.2500 0.3795 0.3000 0.3847 0.3500 0.3909 0.4000 0.3981 0.4500 0.4064 0.5000 0.4158 0.5500 0.4263 0.6000 0.4381 0.6500 0.4511 0.7000 0.4654 0.7500 0.4811 0.8000 0.4982 0.8500 0.5169 0.9000 0.5371 0.9500 0.5590 1.0000 0.5826 图 1精确解2、 用欧拉法求解程序如下：建立函数文件cwfa1.mfunction x,y=cwfa1(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=y(n)+h*feval(fun,x(n),y(n);endx=x'y=y'在MATLAB输入以下程序：>> clear all>> fun=inline(' y*sinx ');>> x,y=cwfa1(fun,0,1,exp(-1),0.05);>> x,y>> plot(x,y,'r+-')ans = 0 0.3679 0.0500 0.3679 0.1000 0.3688 0.1500 0.3706 0.2000 0.3734 0.2500 0.3772 0.3000 0.3819 0.3500 0.3876 0.4000 0.3944 0.4500 0.4023 0.5000 0.4113 0.5500 0.4216 0.6000 0.4332 0.6500 0.4462 0.7000 0.4607 0.7500 0.4768 0.8000 0.4947 0.8500 0.5145 0.9000 0.5364 0.9500 0.56051.0000 0.5871结果及其图象： 图2欧拉解3、用改进欧拉法求解：程序如下：建立函数文件cwfa2.mfunction x,y=cwfa2(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(fun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2;endx=x'y=y'在MATLAB输入以下程序：>> clear all>> fun=inline(' y*sinx'); >> x,y=cwfa2(fun,0,1,exp(-1),0.05);>> x,y>> plot(x,y,'r+-') ans = 0 0.3679 0.0500 0.3683 0.1000 0.3697 0.1500 0.3720 0.2000 0.3753 0.2500 0.3796 0.3000 0.3848 0.3500 0.3911 0.4000 0.3985 0.4500 0.4071 0.5000 0.4169 0.5500 0.4279 0.6000 0.4404 0.6500 0.4544 0.7000 0.4700 0.7500 0.4873 0.8000 0.5066 0.8500 0.5279 0.9000 0.5515 0.9500 0.57761.0000 0.6065结果及其图象： 图3改进欧拉解4、用3阶龙格库塔求解程序如下：建立函数文件cwfa4.mfunction x,y=cwfa4(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=feval(fun,x(n)+h,y(n)+h*(2*k2-k1); y(n+1)=y(n)+h*(k1+4*k2+k3)/6;endx=x'y=y'在Matlab输入以下程序：clear all;fun=inline(' y*sinx');x,y=cwfa4(fun,0,1,exp(-1) ,0.05);x,y plot(x,y, 'b*-')ans = 0 0.3679 0.0500 0.3683 0.1000 0.3697 0.1500 0.3720 0.2000 0.3753 0.2500 0.3796 0.3000 0.3848 0.3500 0.3911 0.4000 0.3985 0.4500 0.4071 0.5000 0.4169 0.5500 0.4280 0.6000 0.4404 0.6500 0.4544 0.7000 0.4700 0.7500 0.4874 0.8000 0.5066 0.8500 0.5280 0.9000 0.5516 0.9500 0.5777 1.0000 0.6065结果及图像： 图4 3阶龙格-库塔解5、 用4阶龙格库塔求解程序如下：建立函数文件cwfa3.mfunction x,y=cwfa3(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); k2=feval(fun,x(n)+h/2,y(n)+h/2*k1);待添加的隐藏文字内容3 k3=feval(fun,x(n)+h/2,y(n)+h/2*k2); k4=feval(fun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x'y=y'在MATLAB输入以下程序：>> clear all;>> fun=inline(' y*sinx ');>> x,y=cwfa3(fun,0,1,exp(-1),0.05);>> x,y >> plot(x,y, 'b*-') ans = 0 0.3679 0.0500 0.3683 0.1000 0.3697 0.1500 0.3720 0.2000 0.3753 0.2500 0.3796 0.3000 0.3848 0.3500 0.3911 0.4000 0.3985 0.4500 0.4071 0.5000 0.4169 0.5500 0.4279 0.6000 0.4404 0.6500 0.4544 0.7000 0.4700 0.7500 0.4874 0.8000 0.5066 0.8500 0.5280 0.9000 0.5516 0.9500 0.5777 1.0000 0.6065结果及其图象： 图5 4阶龙格-库塔解6、精确解与欧拉方法、改进欧拉方法、3阶龙格库塔法以及4阶龙格库塔法得出的解得比较。x=0:0.05:1;y1 = 0.3679 0.3683 0.3697 0.3720 0.3753 0.3795 0.3847 0.3909 0.3981 0.4064 0.4158 0.4263 0.4381 0.4511 0.4654 0.4811 0.4982 0.5169 0.5371 0.5590 0.5826;>> y2= 0.3679 0.3679 0.3688 0.3706 0.3734 0.3772 0.3819 0.3876 0.3944 0.4023 0.4113 0.4216 0.4332 0.4462 0.4607 0.4768 0.4947 0.5145 0.5364 0.5605 0.5871;>> y3= 0.3679 0.3683 0.3697 0.3720 0.3753 0.3796 0.3848 0.3911 0.3985 0.4071 0.4169 0.4279 0.4404 0.4544 0.4700 0.4873 0.5066 0.5279 0.5515 0.5776 0.6065;>> y4= 0.3679 0.3683 0.3697 0.3720 0.3753 0.3796 0.3848 0.3911 0.3985 0.4071 0.4169 0.4280 0.4404 0.4544 0.4700 0.4874 0.5066 0.5280 0.5516 0.5777 0.6065;>>y5= 0.3679 0.3683 0.3697 0.3720 0.3753 0.3796 0.3848 0.3911 0.3985 0.4071 0.4169 0.4279 0.4404 0.4544 0.4700 0.4874 0.5066 0.5280 0.5516 0.5777 0.6065; plot(x,y1,'r+-')hold on,plot(x,y2,'b-'),title('精确解与欧拉方法比较')>> plot(x,y1,'r+-')hold on,plot(x,y3,'b-'),title('精确解与欧拉改进方法比较')>> plot(x,y1,'r+-')hold on,plot(x,y4,'b-'),title('精确解与3阶龙格库塔方法比较')>> plot(x,y1,'r+-')hold on,plot(x,y5,'b-'),title('精确解与4阶龙格库塔方法比较') 图6精确解与欧拉方法比较 图7精确解与欧拉改进方法比较 图8精确解与3阶龙格-库塔方法比较 图9精确解与4阶龙格-库塔方法比较四、结果分析由以上的各种方法与精确解的比较图可以看出在计算精度上，四阶龙格-库塔方法的误差最小，其次是三阶龙格-库塔方法，再次之就是改进欧拉方法，欧拉方法误差则比较大，所以四阶龙格-库塔方法得到最佳的精度。而在计算量上面，相应地，很明显的四阶龙格-库塔方法也是最大，三阶龙格-库塔次之，其后为改进欧拉方法，欧拉方法计算量最小。这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比。我们在实际运用与操作中，可以根据实际情况，选择这4种方法中的其中一种最适合的，追求精度的话，可以使用三阶或四阶龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对的估计上，相应的，精度则有所欠缺。以上的选择，都取决于具体的情况。五、设计心得 这次实验花了较多的时间，先是选择方程方面遇到了困难：有些方程在运行过程中会不断出错，错误的原因有很多，有的是在输入格式上有错，有的是在步长和初值的取法上出错。在试了十多个方程之后终于找到一个自己满意的。接下来的运行中也出现了一些错误。在用第一个方程求精确解时由于没有区分点乘和乘导致错误，有些方程精确解时可以运行但是在用欧拉方程和改进的欧拉方法时就会报错，有的方程即使求的出来结果与精确解的误差是很大的，显然是错误的。还有就是刚开始犯了一个低级错误：把原函数直接代入欧拉求解这样当然是不正确的。取初值对函数图像的影响很大，如果初值取的不合适误差也会非常大，所以要合理的选择初值，在几种求解方法中只有龙格库塔的求解与精确值最为接近几乎是重合的。通过的后面的不同步长求解的比较发现步长越小与精确值越接近。总之在做这次实验中对学习数学有了更多的体会，那就是一定要认真，不能丝毫马虎，一个小小的错误就会导致所有的数据不能运行。六、参考文献1 刘卫国.MATLAB程序设计与应用程. 北京: 高等教育出版社，20092 姜健飞、胡良剑、唐俭.数值分析及其MATLAB实验上海：科学出版社，20043朱功勤 数值计算方法 合肥工业大学出版社。.
