欧拉法解常微分方程

1、实用标准文案大全数学与计算科学学院验报告实验项目名称 Eular方法求解一阶常微分方程数值解所属课程名称偏微分方程数值解实验类型验证性实验日期 2015-3-26班级学号姓名成绩一、实验概述: 【实验目的】熟练掌握应用显性Eular法和隐式Eular法求解一般一阶常微分方程的近似数 值解。【实验原理】虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特 殊类型的方程。求解从实际问题当中归结出来的微分方程主要靠数值解法。欧拉方法是一类重要的数值解法。这类方法回避解y(x)的函数表达式，而是寻求它在一系 列离散节点上的近似值，相邻的两个节点的间距称作步长。假定步长为定数。欧拉方法是

2、一类离散化方法，这类方法将寻求解y(x)的分析问题转化为计算离 散值值的代数问题，从而使问题获得了实质性的简化。然而随之带来的困难是，由于数据量往往很大，差分方法所归结出的可能是个大规模的代数方程组。 【实验环境】1.硬件环境稣处理器安装内存(MM);系统超：苇设低插：系统汾丞询用Intel(R) Core(TM) i3-311OM CPU 240GH? 2.40 GHz4,00 GB (3.85 G3 可用)64位逐f修珠没有可用于此显示器的茎或在控墟人2.2.软件环境Windows 版本Windows 7 揖®!版版院侑 © 2009 Microsoft Corpora

3、tion,保翟所Wt又机Service Pack 1MATLAB7.0二、实验内容：【实验过程】(实验步骤)(一)实验任务描述某种化学反应过程的方程，利用显性和隐形Eualar方法求解下列一阶线性微分方程组的近似数值解：d%=-0.04yi +10丫2jdy%=0.04yi+l04yiy2-3Ml07y2dy =3-y2 dt.(0)=1, y2(0)=Q y3(0)=0(二)求解过程Eular方法：一阶线性微分方程初值问题7= f (x,y),a <x <b <、y(a) = y0a=xo<x1<.<xn=b (1)Xn = Xo nh, h为步长方程离散

4、化：差分和差商yi -yoyi - yOy (xo)：Xi -xohf(xo, yo) =* J0 h yi =y0 +hf (x°,yo)(2)Vn 1 =yn hf(x0,yo)通过初始值yo,依据递推公式(2)逐步算出yi,y2,., yn就为显性的Eular方 法。隐形Eular方法:（3）、=y0 +hf（x1，yi）=Jn4t=yn +hf Un由，Yn 十）公式(3)即为隐式Eular公式。(三)程序算法1.利用显式Eular法方求解利用MATLA进行求解，编写脚本文件如下:文件名：hql.m% 显性Eular方法f0=1; g0 =0;z0=0delta=0.01;t

5、ime=1;t=0:delta:time;f=zeros(size(t);g=zeros(size(t);z=zeros(size(t);f1=zeros(size(t);g1=zeros(size(t);z1=zeros(size(t);f(1)=f0;g(1)=g0；z(1)=z0;for i=2:length(t)f1(i-1) = -0.04*f(i-1) + 10000*f(i-1)*g(i-1);f(i)=f(i-1)+f1(i-1)*delta;g1(i-1) = 0.04*f(i-1) - 10000*f(i-1)*g(i-1)-3*10A7*g(i-1)A2;g(i)=g(i

6、-1)+g1(i-1)*delta；z1(i-1)=3*10A7*g(i-1)A2;z(i)=z(i-1)+z(i-1)*delta;Fun=f+g+zendfigureplot(t,f, 'o')；xlabel('t');ylabel('y1');title( 't-y1 变化图')figureplot(t,g, 'o');xlabel('t');ylabel('y2');title( 't-y2 变化图')figureplot(t,z, 'o');

7、xlabel('t');ylabel('y3');title( 't-y3 变化图')figureplot(t,Fun);xlabel('t');ylabel('y1+y2+y3');title( 't-y1+y2+y3 变化图)【实验结论】A步长h=0.001时进行数据测试。结果如下:迭代第一次时,l.OOQOe-* 1x101 de0Columns 12 through 22结果与方程描述内容相符。迭代第二次时,Columns I through 111.0000 L 0000E00000.9959结果

8、与方程描述内容基本相符。迭代三次时,Fun -Columns 1 thraugh 11EOOOO k 0000 L 00000.99590.980800八-alColumns 12 through 22结果与方程描述内容基本相符。5-Q5(H)Q也MM也MM-0, D0007, WOO-0.0M03MlNaXFlaNNaNMMTJaNNaffNaJiNaffNaifMM迭代1000次时,Fun =和 3 *Caliums I through i I6 MM D.g 钟 MMgColumns 12 thraufh 22二工 lfi-55 -InfColiuicas 23 through 33Hm

9、兵 Colujufts 3l through 44模拟结果已经严重脱离事实，故当选择delta为0.001时，该迭代方法不收敛。时间与个变量直接的变化关系如图所示:从上述图形可以明显看出，在迭代的不断进行时，各变量与时间的变化越来越 大，且严重脱离了方程所描述的现实意义。B.当选择h=0.00000001时，模拟结果如下:迭代第一次,Fun =Columns 1 through IE1100Columns 19 through 350000与A中结果相同。Fun =Colwhns1Columns迭代第二次，1 through 1811019 through 36跌二次迭代结果明显优于一中。跌三

10、次迭代结果,Fun -Columns值l.OOOOe-A- 1x101 di1"01 de1x1 1 10JColumns1 through 1811 I19 through 36并未产生误差。地1000次迭代结果,Colums 1 thxouEh 10I.OWl.omLQOQQLOWL MML WQL QMQLMML MQQColuw 11 thrcniEh 20lmmLornl omLornLmolmqqlmq。i , 0000LQggCdliwH 21 thrwili JO构AuomLom Lmolmqq lmw l qoooColiwis. 31 thmugh 40lowI,

11、l GmLornk«w。lwoolmel, coool oooo1shihdt W1 throuih 5Ql.omIRDWJLQQMLDmIL WMQLMQQL(NMI, MMLQQgL QEQCglvims 5：l tkiuEh <501通1.八fiJOJ dJ |311有巾J泰1 mi Jpl 3A dt a 上de V >ft结果明显是收敛的。时间与个变量直接的变化关系如图所示:t-y1 +y2+y凌化图00.10.20.30.40.50,60 70.80.91t乂d从图中能够清晰看出，当 h=0.00000001时，模拟结果与方程所表示的显示意义相吻合。说明了显性

12、Eualr方法的收敛性是与步长的选择是相关。这就对我们们选择步长造成了困难， 由于选择的步长不合适有可能得出错误的结论。【实验小结】（收获体会）1、软件使用在写MATLABS言的时候要深刻理解题的意图，整理好思绪再做题目，在我运算 的过程中，h取值取得越小、越细微，曲线逼近的越好。2、欧拉法的缺点简单地取切线的端点作为下一的起点进行计算，当步数增多时，误差会因积累 而越来越大。因此欧拉格式一般不用于实际计算。3、实验感想在这次上机实验中，我掌握了解决常微分方程的基本方法，同时学会使用计算机 软件对两种不同方法得到的结果进行判断，对我们以后对数据进行分析很有祁助。三、指导教师评语及成绩：评语评语

13、等级优良中及 格/、及格1.实验报告按时完成，字迹清楚，文字叙述流畅，逻辑性强2.实验方案设计合理3.实验过程（实验步骤详细，记录完整，数据合理，分析透彻）4实验结论正确.成绩：指导教师签名： 批阅日期：附录：源程序程序1:%&Ti Eular 方法f0=1; g0 =0;z0=0delta=0.00000001;time=0.00001;t=0:delta:time;f=zeros(size(t);g=zeros(size(t);z=zeros(size(t);f1=zeros(size(t);g1=zeros(size(t);z1=zeros(size(t);f(1)=f0;g(1

14、)=g0；z(1)=z0;for i=2:length(t)f1(i-1) = -0.04*f(i-1) + 10000*f(i-1)*g(i-1);f(i)=f(i-1)+f1(i-1)*delta;g1(i-1) = 0.04*f(i-1) - 10000*f(i-1)*g(i-1)-3*10A7*g(i-1)A2;g(i)=g(i-1)+g1(i-1)*delta；z1(i-1)=3*10A7*g(i-1)A2;z(i)=z(i-1)+z(i-1)*delta;Fun=f+g+zendfigureplot(t,f, 'o')；xlabel( 't');ylabel( 'y1');title( 't-y1 变化图)figureplot(t,g, 'o')；xlabel( 't');ylabel( '