资源描述:
《一维有限差分法之稳定性分析》由会员上传分享,免费在线阅读,更多相关内容在工程资料-天天文库。
1、2013~2014学年第一学期实验报告班级 学号:姓名: 实验时间:11月29日 成绩:实验项目一维有限差分法之稳定性分析所属课程微分方程数值解实验目的了解一维有限差分法之算法构造,如Euler法,改进Euler等。实验环境Matlab7.0实验内容掌握一维有限差分法之算法构造,用Euler法、Taylor法求解课本上45页实习题,并分析稳定性实验过程第二题:(1)利用Euler法计算出U(2),用Milne迭代法+Jacobi法:定义函数f2functiony=f2(t,u)y=t^3*u^3-t*u;%Milne迭代法+Jacobi方法,例2%Mil
2、ne迭代法+Jacobi迭代法,例2n=input('n=');%n等分T=input('T=');%时间上界H=input('H=');%小区间H等分IT=input('IT=');%Jacobi迭代次数,内迭代次数h=T/n;hh=h/H;u=zeros(n+1,1);u(1)=0.5;x(1)=0.5;formm=1:Hx(mm+1)=x(mm)+hh*f2(mm*hh,x(mm));endu(2)=x(H+1);form=1:n-1%Milne开始u0=u(m)+h*f2(m*h,u(m));fork=1:IT%Jacobi%Jacobi迭代开始
3、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结束%u'7t=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*-')holdonu2=zeros(n+1,1);u2(1)=0.5;x(1)=0.5;formm=1:Hx(mm+1)=x(mm)+hh
4、*f2(mm*hh,x(mm));endu2(2)=x(H+1);form=1:n-1%Milne开始u0=u2(m)+h*f2(m*h,u2(m));fork=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*h^3*u0^2-(m+2)*h));u0=u1;end%Jacobi%Jacobi迭代结束u2(m+2)=u0;end%Milne结束%u't=0:
5、h:T;uu=1./sqrt(t.^2+1+3*exp(t.^2));subplot(2,1,2)u2'-uuplot(t,u2,'go-',t,uu,'r*-')结果:(画图如下)(2)用Adams二步内插法+Newton法:n=input('n=');T=input('T=');H=input('H=');7IT=input('IT=');h=T/n;hh=h/H;u=zeros(n+1,1);u(1)=0.5;x(1)=0.5;formm=1:Hx(mm+1)=x(mm)+hh*f2(mm*hh,x(mm));endu(2)=x(H+1);form
6、=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));fork=1:ITu1=u0-(u0-5/12*h*f2((m+2)*h,u0)+w)/(1-5/12*h*(3*(m+2)^3*h^3*u0^2-(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>-')结果:第三题:定义函数f3functiony=f
7、3(x,y)y=x-y-x*(x^2+y^2);functiony=g3(x,y)y=x+y-y*(x^2+y^2);%Milne迭代法+Jacobi方法,例3n=input('n=');%n等分7T=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;y1(1)=1;formm=1:Hx1(mm+1)=x1(mm)+hh*f3(mm*h
8、h,x1(mm));y1(mm+1)=y1(mm)+hh*g3(mm*hh,y1