实验四 常微分方程初值问题数值解法.doc

佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 常微分方程初值问题数值解法 专业班级 12.数学与应用数学（师范） 姓名 叶楚欣 学号 2012214103 指导教师 黄国顺 成 绩 日 期 一. 实验目的1、 理解如何在计算机上实现用Euler法、改进Euler法、RungeKutta算法求一阶常微分方程初值问题 的数值解。2、 利用图形直观分析近似解和准确解之间的误差。3、 学会Matlab提供的ode45函数求解微分方程初值问题。二、实验要求（1） 按照题目要求完成实验内容；（2） 写出相应的Matlab 程序；（3） 给出实验结果(可以用表格展示实验结果)；（4） 分析和讨论实验结果并提出可能的优化实验。（5） 写出实验报告。三、实验步骤1、用编好的Euler法、改进Euler法计算书本P167 的例1、P171例题3。（1）取，求解初值问题（2）取，求解初值问题2、用RungeKutta算法计算P178例题、P285实验任务（2）（1）取，求解初值问题（2）求初值问题的解在处的近似值，并与问题的解析解相比较。3、用Matlab绘图函数plot(x,y)绘制P285实验任务（2）的精确解和近似解的图形。4、使用matlab中的ode45求解P285实验任务（2），并绘图。四、实验结果1、Euler算法程序、改进Euler算法程序；2、用Euler算法程序、改进Euler算法求解P167例题1的运行结果；3、RungeKutta算法程序；4、用RungeKutta算法求解P178例题、P285实验任务（2），计算结果如下（其中表示数值解，表示解析解，结果保留八位有效数字）：0.050.10.150.20.25-0.02219009-0.03877057-0.04975651-0.05516258-0.05500309-0.02219009-0.03877058-0.04975651-0.05516258-0.055003100.30.350.40.450.5-0.04929202 -0.03804297-0.021269240.001016230.02880079-0.04929202-0.03804298-0.021269250.001016220.028800785、P285实验任务（2）精确解与近似解的图形比较（图片贴到此处）6、用matlab中的ode45求解P285实验任务（2） MATLAB desktop keyboard shortcuts, such as Ctrl+S, are now customizable. In addition, many keyboard shortcuts have changed for improved consistency across the desktop. To customize keyboard shortcuts, use Preferences. From there, you can also restore previous default settings by selecting R2009a Windows Default Set from the active settings drop-down list. For more information, see Help. Click here if you do not want to see this message again. ydot_fun=inline(y-2*x./y,x,y); x,y=euler_f(ydot_fun,0,1,0.1,10)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 8 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 Columns 9 through 11 0.800000000000000 0.900000000000000 1.000000000000000y = Columns 1 through 4 1.000000000000000 1.100000000000000 1.191818181818182 1.277437833714722 Columns 5 through 8 1.358212599560289 1.435132918657796 1.508966253566332 1.580338237655217 Columns 9 through 11 1.649783431047711 1.717779347860087 1.784770832497982 ydot_fun=inline(y-2*x./y,x,y); x,y=euler_r(ydot_fun,0,1,0.1,10)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 8 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 Columns 9 through 11 0.800000000000000 0.900000000000000 1.000000000000000y = Columns 1 through 4 1.000000000000000 1.095909090909091 1.184096569242997 1.266201360875776 Columns 5 through 8 1.343360151483999 1.416401928536909 1.485955602415669 1.552514091326146 Columns 9 through 11 1.616474782752058 1.678166363675186 1.737867401035414 ydot_fun=inline(-y+x+1,x,y); x,y=euler_f(ydot_fun,0,1,0.1,5)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 6 0.400000000000000 0.500000000000000y = Columns 1 through 4 1.000000000000000 1.000000000000000 1.010000000000000 1.029000000000000 Columns 5 through 6 1.056100000000000 1.090490000000000 ydot_fun=inline(-y+x+1,x,y); x,y=euler_r(ydot_fun,0,1,0.1,5)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 6 0.400000000000000 0.500000000000000y = Columns 1 through 4 1.000000000000000 1.005000000000000 1.019025000000000 1.041217625000000 Columns 5 through 6 1.070801950625000 1.107075765315625 ydot_fun=inline(y2,x,y); x,y=Runge_Kutta4(ydot_fun,0,1,0.1,5)x = Columns 1 through 4 0 0.100000000000000 0.200000000000000 0.300000000000000 Columns 5 through 6 0.400000000000000 0.500000000000000y = Columns 1 through 4 1.000000000000000 1.111110490052194 1.249997992047015 1.428566186301445 Columns 5 through 6 1.666653257250323 1.999963258950669 ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10)x = Columns 1 through 4 0 0.050000000000000 0.100000000000000 0.150000000000000 Columns 5 through 8 0.200000000000000 0.250000000000000 0.300000000000000 0.350000000000000 Columns 9 through 11 0.400000000000000 0.450000000000000 0.500000000000000y = Columns 1 through 4 0 -0.022190087076823 -0.038770573733692 -0.049756511058554 Columns 5 through 8 -0.055162578526671 -0.055003093175772 -0.049292018554665 -0.038042973450910 Columns 9 through 11 -0.021269240403008 0.001016225997586 0.028800791009418ydot_fun=inline(-y+x2+4*x-1)./2,x,y);x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10); sprintf(%.8f,%.8fn,x,y)ans =0.00000000,0.050000000.10000000,0.150000000.20000000,0.250000000.30000000,0.350000000.40000000,0.450000000.50000000,0.00000000-0.02219009,-0.03877057-0.04975651,-0.05516258-0.05500309,-0.04929202-0.03804297,-0.021269240.00101623,0.02880079xx=0:0.05:0.5;z=exp(-xx./2)+xx.2-1;sprintf(%.8f,%.8fn,xx,z)ans =0.00000000,0.050000000.10000000,0.150000000.20000000,0.250000000.30000000,0.350000000.40000000,0.450000000.50000000,0.00000000-0.02219009,-0.03877058-0.04975651,-0.05516258-0.05500310,-0.04929202-0.03804298,-0.021269250.00101622,0.02880078 xx=0:0.05:0.5;z=exp(-xx./2)+xx.2-1;hold on plot(x,y,o);plot(xx,z,*); hold off ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=ode45(ydot_fun,0,0.5,0); x;yans = Columns 1 through 4 0 0.000100475457260 0.000200950914521 0.000301426371781 0 -0.000050226371419 -0.000100430028501 -0.000150610971371 Columns 5 through 8 0.000401901829042 0.000904279115343 0.001406656401645 0.001909033687947 -0.000200769200158 -0.000451219637267 -0.000701102241287 -0.000950417028058 Columns 9 through 12 0.002411410974249 0.004923297405759 0.007435183837268 0.009947070268778 -0.001199164013417 -0.002434382472993 -0.003655408270323 -0.004862243380400 Columns 13 through 16 0.012458956700288 0.024958956700288 0.037458956700288 0.049958956700288 -0.006054889775763 -0.011778983079828 -0.017151998201403 -0.022174175468426 Columns 17 through 20 0.062458956700288 0.074958956700288 0.087458956700288 0.099958956700288 -0.026845753781789 -0.031166970574532 -0.035138061683373 -0.038759261501758 Columns 21 through 24 0.112458956700288 0.124958956700288 0.137458956700288 0.149958956700288 -0.042030803032876 -0.044952917848809 -0.047525835962155 -0.049749785978392 Columns 25 through 28 0.162458956700288 0.174958956700288 0.187458956700288 0.199958956700288 -0.051624995148617 -0.053151689328665 -0.054330092850809 -0.055160428675467 Columns 29 through 32 0.212458956700288 0.224958956700288 0.237458956700288 0.249958956700288 -0.055642918443663 -0.055777782436119 -0.055565239445031 -0.055005506925134 Columns 33 through 36 0.262458956700288 0.274958956700288 0.287458956700288 0.299958956700287 -0.054098801045891 -0.052845336650564 -0.051245327128055 -0.049298984563352 Columns 37 through 40 0.312458956700288 0.324958956700288 0.337458956700287 0.349958956700287 -0.047006519789452 -0.044368142346404 -0.041384060353229 -0.038054480657744 Columns 41 through 44 0.362458956700287 0.374958956700288 0.387458956700288 0.399958956700287 -0.034379608888238 -0.030359649412486 -0.025994805209740 -0.021285278019953 Columns 45 through 48 0.412458956700287 0.424958956700287 0.437458956700287 0.449958956700287 -0.016231268395211 -0.010832975658719 -0.005090597776859 0.000995668492165 Columns 49 through 52 0.462458956700287 0.471844217525216 0.481229478350144 0.490614739175072 0.007425627547137 0.012479158933356 0.017726249535970 0.023166817935270 Column 53 0.500000000000000 0.028800783076419 plot(x,y,*)五、讨论分析（通过最后绘制的图形直观分析近似解与准确解之间的误差，自己补充）从图中可以看到4阶的Runge-Kutta画出来的点基本和原函数的对应点相重合，精度之高，从数据来看，起码有9位有效数字。六、改进实验建议（自己补充）可以通过输出数据的精度控制函数实验四 常微分方程初值问题数值解法（这里只提供解答格式请同学自己按照上机实验报告格式填写）1、饿uler法求解初值问题function x,y=euler_f(ydot_fun, x0, y0, h, N)% Euler（向前）公式，其中% ydot_fun - 一阶微分方程的函数% x0, y0 - 初始条件% h - 区间步长% N - 区间的个数% x - Xn 构成的向量% y - Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(ydot_fun, x(n), y(n);end用欧拉法计算p167例1ydot_fun=inline(y-2*x./y,x,y); x,y=euler_f(ydot_fun,0,1,0.1,10)x = Columns 1 through 110 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.5000000000000000.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y =.000000000000000 1.100000000000000 1.191818181818182 1.277437833714722 1.358212599560289 1.4351329186577961.508966253566332 1.580338237655217 1.649783431047711 1.717779347860087 1.7847708324979822、改进Euler公式function x,y=euler_r(ydot_fun, x0, y0, h, N)% 改进Euler公式，其中% ydot_fun - 一阶微分方程的函数% x0, y0 - 初始条件% h - 区间步长% N - 区间的个数% x - Xn 构成的向量% y - Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot_fun, x(n), y(n); y(n+1)=y(n)+h/2*(feval(ydot_fun, x(n), y(n)+feval(ydot_fun, x(n+1), ybar);end用改进欧拉公式计算P167例题1ydot_fun=inline(y-2*x./y,x,y); x,y=euler_r(ydot_fun,0,1,0.1,10)x =Columns 1 through 60 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000Columns 7 through 110.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.000000000000000y =Columns 1 through 61.000000000000000 1.095909090909091 1.184096569242997 1.266201360875776 1.343360151483999 1.416401928536909Columns 7 through 111.485955602415669 1.552514091326146 1.616474782752058 1.678166363675186 1.737867401035414用欧拉法、改进欧拉法计算P171例题3ydot_fun=inline(-y+x+1,x,y); x,y=euler_f(ydot_fun,0,1,0.1,5)x =0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000y =1.000000000000000 1.000000000000000 1.010000000000000 1.029000000000000 1.056100000000000 1.090490000000000用改进欧拉法计算P171例题3ydot_fun=inline(-y+x+1,x,y); x,y=euler_r(ydot_fun,0,1,0.1,5)x=0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000y=1.000000000000000 1.005000000000000 1.019025000000000 1.041217625000000 1.070801950625000 1.1070757653156253、标准四阶Runge_Kutta公式function x,y=Runge_Kutta4(ydot_fun, x0, y0, h, N)% 标准四阶Runge_Kutta公式，其中% ydot_fun - 一阶微分方程的函数% x0, y0 - 初始条件% h - 区间步长% N - 区间的个数% x - Xn 构成的向量% y - Yn 构成的向量x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0;for n=1:N x(n+1)=x(n)+h; k1=h*feval(ydot_fun, x(n), y(n); k2=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k1); k3=h*feval(ydot_fun, x(n)+1/2*h, y(n)+1/2*k2); k4=h*feval(ydot_fun, x(n)+h, y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end用四阶龙格库塔计算P178例题 ydot_fun=inline(y2,x,y); x,y=Runge_Kutta4(ydot_fun,0,1,0.1,5)x =0 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000y =1.000000000000000 1.111110490052194 1.249997992047015 1.428566186301445 1.666653257250323 1.999963258950669用龙格库塔方法计算P285实验任务（2）ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=Runge_Kutta4(ydot_fun,0,0,0.05,10)x =Columns 1 through 60 0.050000000000000 0.100000000000000 0.150000000000000 0.200000000000000 0.250000000000000Columns 7 through 110.300000000000000 0.350000000000000 0.400000000000000 0.450000000000000 0.500000000000000y = Columns 1 through 60 -0.022190087076823 -0.038770573733692 -0.049756511058554 -0.055162578526671 -0.055003093175772Columns 7 through 11-0.049292018554665 -0.038042973450910 -0.021269240403008 0.001016225997586 0.028800791009418xx=0:0.05:0.5; z=exp(-xx./2)+xx.2-1; hold on plot(x,y,o);plot(xx,z,*)hold off其图形如图一所示图一 精确解与龙格库塔计算结果比较4、 用Matlab提供的ode45求解P285实验任务（2）ydot_fun=inline(-y+x2+4*x-1)./2,x,y); x,y=ode45(ydot_fun, 0,0.5, 0); x;yplot(x,y,*)ans = Columns 1 through 60 0.000100475457260 0.000200950914521 0.000301426371781 0.000401901829042 0.000904279115343 0 -0.000050226371419 -0.000100430028501 -0.000150610971371 -0.000200769200158 -0.000451219637267 0.001406656401645 0.001909033687947 0.002411410974249 0.004923297405759 0.007435183837268 0.009947070268778 -0.000701102241287 -0.000950417028058 -0.001199164013417 -0.002434382472993 -0.003655408270323 -0.0048622433804000.012458956700288 0.024958956700288 0.037458956700288 0.049958956700288 0.062458956700288 0.074958956700288 -0.006054889775763 -0.011778983079828 -0.017151998201403 -0.022174175468426 -0.026845753781789 -0.03116