实验一常微分方程的数值解法.docxVIP

《偏微分方程数值解》实验报告实验一常微分方程的数值解法实验内容分别用欧拉法和改进的欧拉法（预报—校正格式）和四阶龙格库塔法求解初值问题。（h=0.1） (1)并与真解作比较（列出表格，并画图）选取一种理论上收敛但是不稳定的算法对问题（1）进行计算，并与真解作比较（列表，画图）。实验结果分析1、欧拉法第一步：编写Euler.m函数，代码如下function Euler(fun,x0,y0,xn,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xn-x0)/N;for n=1:Nx(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n));endT=[x,y]y1 = dsolve(Dy=-y+x+1,y(0)=1,x)x1 = 0:0.1:1;y2=x1 + exp(-x1);A=[x1,y2]plot(x,y,O,x1,y2);xlabel(x);ylabel(y(x));title(??êμ?aó??·à-·¨?ó?a??±è);gridon;legend(?·à-?a,??êμ?a);Euler(fun,0,1,1,10)其中fun.m为function z=fun(x,y)z=-y+x+1;程序结果：2、改进的欧拉法编写TranEuler.m函数，代码如下functionTranEuler(fun,x0,y0,xn,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xn-x0)/N;for n=1:Nx(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n)); y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));endT=[x,y]y1 = dsolve(Dy=-y+x+1,y(0)=1,x)x1 = 0:0.1:1;y2=x1 + exp(-x1);A=[x1,y2]plot(x,y,O,x1,y2);xlabel(x);ylabel(y(x));title(??êμ?aó?????μ??·à-·¨?ó?a??±è);gridon;legend(?????·à-?a,??êμ?a);程序结果：3、4阶Rk法编写4阶Rk函数，代码如下clearall; clc;ticfun=inline(-y+x+1,x,y);[x,y]=ode45(fun,[0,1],1);B=[x,y]y1 = dsolve(Dy=-y+x+1,y(0)=1,x)x1 = 0:0.1:1;y2=x1 + exp(-x1);A=[x1,y2]plot(x,y,O,x1,y2);xlabel(x);ylabel(y(x));title(??êμ?aó?4?×RK·¨?ó?a??±è);gridon;legend(4?×RK·¨?a,??êμ?a);toc程序结果：010.02500000000000001.000309908082360.05000000000000001.001229418796880.07500000000000001.002743483983400.1000000000000001.004837418333330.1250000000000001.007496899304160.1500000000000001.010707971546880.1750000000000001.014457018923110.2000000000000001.018730753616130.2250000000000001.023516216053550.2500000000000001.028800778913400.2750000000000001.034572121804180.3000000000000001.040818221412120.3250000000000001.047527351431190.3500000000000001.054688086187990.3750000000000001.062289277731270.4000000000000001.070320046916840.4250000000000001.078769783344220.4500000000000001.087628148636600.4750000000000001.096885055710540.5000000000000001.106530660709310.5250000000000001.116555362945530.5500000000000001.126949807869