一维有限差分法之稳定性分析报告

《一维有限差分法之稳定性分析报告》由会员分享，可在线阅读，更多相关《一维有限差分法之稳定性分析报告（12页珍藏版）》请在装配图网上搜索。

1、范文 范例指导 参考20 1320 14学年第一 学期实验报告班级学号：姓名：实验时间： 11 月 29 日成绩：实验一维有限差分法之稳定性分析所属微分方程数值解项目实验目的实验环境实验性内容课程了解一维有限差分法之算法构造，如 Euler 法，改进 Euler 等 。Matlab 7.0掌握一维有限差分法之算法构造，用 Euler 法、 Taylor 法求解课本上45 页实习题 ，并分析稳定第二题 ：（ 1）利用 Euler 法计算出U（ 2）， 用 Milne 迭代法 +Jacobi 法：实定义函数f2验 function y=f2(t,u) y=t3*u3-t*u;过 % Milne 迭

2、代法 +Jacobi 方法，例 2%Milne迭代法 +Jacobi 迭代法 ，例 2程T=input(T=);%时间上界H=input(H=);%小区间 H 等分word 资料整理分享范文 范例指导 参考IT=input(IT=);%Jacobi迭代次数 ，内迭代次数h=T/n;hh=h/H;u=zeros(n+1,1);u(1)=0.5;x(1)=0.5;for mm=1:Hx(mm+1)=x(mm)+hh*f2(mm*hh,x(mm);endu(2)=x(H+1);for m=1:n-1%Milne开始u0=u(m)+h*f2(m*h,u(m);for k=1:IT%Jacobi%Jac

3、obi 迭代开始w=-u(m)-4/3*h*f2(m+1)*h,u(m+1)-1/3*h*f2(m*h,u(m);u1=1/3*h*f2(m+2)*h,u0)-w;u0=u1;end%Jacobi %Jacobi 迭代结束u(m+2)=u0;end%Milne结束%ut=0:h:T;uu=1./sqrt(t.2+1+3*exp(t.2);subplot(2,1,1)u-uuplot(t,u,bo-,t,uu,r*-)hold onu2=zeros(n+1,1);u2(1)=0.5;x(1)=0.5;for mm=1:Hx(mm+1)=x(mm)+hh*f2(mm*hh,x(mm);endu2(

4、2)=x(H+1);for m=1:n-1%Milne开始u0=u2(m)+h*f2(m*h,u2(m);for k=1:IT %Jacobi 迭代开始w=-u2(m)-4/3*h*f2(m+1)*h,u2(m+1)-1/3*h*f2(m*h,u2(m);u1=u0-(u0-1/3*h*f2(m+2)*h,u0)+w)/(1-1/3*h*(3*(m+2)3*h3*u02-(m+2)*h);u0=u1;word 资料整理分享范文 范例指导 参考end%Jacobi %Jacobi 迭代结束u2(m+2)=u0;end%Milne结束%ut=0:h:T;uu=1./sqrt(t.2+1+3*exp

5、(t.2);subplot(2,1,2)u2-uuplot(t,u2,go-,t,uu,r*-)结果 ：（画图如下 ）0.50.450.40.350.30.2500.10.20.30.40.50.60.70.80.91（ 2）用 Adams 二步内插法 +Newton 法：n=input( n= );T=input( T= );H=input( H= );IT=input( IT= );h=T/n;hh=h/H;u=zeros(n+1,1);word 资料整理分享范文 范例指导 参考u(1)=0.5;x(1)=0.5;for mm=1:Hx(mm+1)=x(mm)+hh*f2(mm*hh,x(

6、mm);endu(2)=x(H+1);for m=1:n-1u0=u(m)+h*f2(m*h,u(m);w=-u(m+1)-2/3*h*f2(m+1)*h,u(m+1)-1/12*h*f2(m*h,u(m);for k=1:ITu1=u0-(u0-5/12*h*f2(m+2)*h,u0)+w)/(1-5/12*h*(3*(m+2)3*h3*u02-(m+2)*h);u0=u1;endu(m+2)=u0;endt=0:h:T;uu=1./sqrt(t.2+1+3*exp(t.2);plot(t,u, r*- ,t,uu, b- )结果 ：word 资料整理分享范文 范例指导 参考0.50.450

7、.40.350.30.2500.10.20.30.40.50.60.70.80.91第三题 ：定义函数 f3function y=f3(x,y)y=x-y-x*(x2+y2);function y=g3(x,y)y=x+y-y*(x2+y2);% Milne迭代法 +Jacobi 方法，例 3n=input(n=);%n等分T=input(T=);%时间上界H=input(H=);%小区间 H 等分IT=input(IT=);%Jacobi迭代次数 ，内迭代次数h=T/n;hh=h/H;x=zeros(n+1,1);y=zeros(n+1,1);x(1)=2;y(1)=1;x1(1)=2;wo

8、rd 资料整理分享范文 范例指导 参考y1(1)=1;for mm=1:Hx1(mm+1)=x1(mm)+hh*f3(mm*hh,x1(mm);y1(mm+1)=y1(mm)+hh*g3(mm*hh,y1(mm);end x(2)=x1(H+1);y(2)=y1(H+1);for m=1:n-1x0=x(m)+h*f3(m*h,x(m);w1=-x(m)-4/3*h*f3(m+1)*h,x(m+1)-1/3*h*f3(m*h,x(m);y0=y(m)+h*g3(m*h,y(m);w2=-y(m)-4/3*h*g3(m+1)*h,y(m+1)-1/3*h*g3(m*h,y(m);for k=1:

9、IT %Jacobi 迭代开始x1=h/3*f3(m+2)*h,x0)-w1;y1=h/3*g3(m+2)*h,y0)-w2;x0=x1;y0=y1;end %Jacobi 迭代结束x(m+2)=x0;y(m+2)=y0;endxy结果：word 资料整理分享范文 范例指导 参考n=100T=1H=10IT=100ans =2.00001.97991.95921.93901.91831.89801.87721.85691.83611.81591.79521.77511.75451.73461.71411.69441.67421.65471.63471.61551.59581.57691.557

10、51.53891.51991.50171.48301.46511.44691.42941.41151.39441.37691.36021.34311.32691.31021.29431.27811.26271.24681.23181.21631.20171.18671.17251.15781.14411.12981.11641.10251.08951.07601.06331.05011.03791.02501.01311.00050.98890.97670.96540.95350.94250.93080.92010.90870.89820.88710.87690.86590.85590.845

11、20.83540.82490.81530.80490.79540.78520.77590.76580.75670.74670.73760.72770.71880.70900.70010.69030.68150.67170.66300.65320.64450.63470.62600.61620.60740.59760.58880.5789ans =word 资料整理分享范文 范例指导 参考1.00001.00011.00041.00061.00111.00151.00221.00271.00361.00421.00521.00601.00701.00801.00911.01011.01141.0

12、1251.01381.01501.01641.01761.01911.02041.02191.02321.02471.02611.02761.02901.03061.03191.03351.03491.03641.03781.03941.04071.04221.04351.04501.04631.04771.04891.05041.05151.05291.05401.05531.05631.05761.05851.05971.06061.06171.06241.06351.06411.06511.06571.06651.06701.06771.06811.06881.06911.06961.0

13、6981.07021.07031.07061.07061.07081.07061.07071.07041.07041.07001.06991.06941.06911.06851.06811.06741.06691.06601.06541.06441.06371.06261.06171.06051.05951.05821.05711.05561.05441.05281.05151.04981.0484（ 2） %Milne 方法 +Newton 法n=input(n=);%n等分T=input(T=);%时间上界H=input(H=);%小区间 H 等分IT=input(IT=);% Newto

14、n迭代次数 ，内迭代次数h=T/n;hh=h/H;x=zeros(n+1,1);word 资料整理分享范文 范例指导 参考y=zeros(n+1,1);x(1)=2;y(1)=1;x1(1)=2;y1(1)=1;for mm=1:H%Euler开始x1(mm+1)=x1(mm)+hh*f3(mm*hh,x1(mm);y1(mm+1)=y1(mm)+hh*g3(mm*hh,y1(mm);end%Euler 结束x(2)=x1(H+1);y(2)=y1(H+1);for m=1:n-1x0=x(m)+h*f3(m*h,x(m);w=-x(m)-4/3*h*f3(m+1)*h,x(m+1)-1/3*

15、h*f3(m*h,x(m);y0=y(m)+h*g3(m*h,y(m);w1=-y(m)-4/3*h*g3(m+1)*h,y(m+1)-1/3*h*g3(m*h,y(m);for k=1:IT%Newton开始x1=x0-(x0-1/3*h*f3(m+2)*h,x0)+w1)/(1-1/3*h*(3*(m+2)3*h3*x02-(m+2)*h);y1=y0-(y0-1/3*h*g3(m+2)*h,y0)+w2)/(1-1/3*h*(3*(m+2)3*h3*y02-(m+2)*h);x0=x1;y0=y1;end%Newton结束x(m+2)=x0;word 资料整理分享实验结果及分析范文 范例

16、指导 参考y(m+2)=y0;endxy结果 :n=100T=1H=10IT=100ans =2.00001.97990.99700.99581.04401.04421.04441.04461.04471.04491.04501.04521.04531.04551.04561.04581.04591.04601.04611.04621.04641.04651.04661.04671.04681.04681.04691.04701.04711.04711.04721.04731.04731.04741.04741.04751.04751.04751.04751.04761.04761.04761

17、.04761.04761.04761.04761.04761.04761.04751.04751.04751.04741.04741.04731.04731.04721.04711.04711.04701.04691.04681.04671.04661.04651.04641.04631.04621.04611.04591.04581.04571.0455word 资料整理分享范文范例指导参考1.04541.04521.04501.04491.04471.04451.04431.04411.04391.04371.04351.04331.04311.04291.04261.04241.0422

18、1.04191.04171.04141.04111.04091.04061.04031.04001.03971.03941.03911.0388ans =1.00001.00011.04871.04871.04871.04881.04881.04881.04891.04891.04891.04891.04901.04901.04901.04901.04901.04911.04911.04911.04911.04911.04921.04921.04921.04921.04921.04921.04931.04931.04931.04931.04931.04931.04931.04931.04931

19、.04931.04931.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04941.04931.04931.04931.04931.04931.04931.04931.04931.04931.04931.04931.04921.04921.04921.04921.04921.04921.04911.04911.04911.04911.04911.04911.04901.04901.04901.04901.04891.04891.04891.04891.04881.04881.04881.04881.04871.04871.04871.04861.04861.04861.04851.04851.04851.04841.0484word 资料整理分享范文 范例指导 参考实通过本次实验掌握一维有限差分法之算法构造，用 Euler 法、 Milne法等。了解了当 n 越大 ， H 越验大时 ，也就是 h 趋于 0 时，误差都会变得更小，稳定性也更好。心得总结教师评语注：以上各栏若填写不够 ，可自行扩展 。word 资料整理分享