姓名_学号_中国海洋大学偏微分方程的数值解法第四次作业

1、偏微分方程的数值解法上机习题四题目：一、PMECME算法1）用四步四阶Adams方法预估（课本P14,k=3）,用三步四阶隐式Adams方法校正（课本P16,k=3）,编写PMECME算法。2）用四步四阶Adams方法预估（课本P14,k=3）,用三步四阶Milne方法校正（课本P17,（1.2.31）,编写PMECME算法。二、变步计算1）用误差的后验估计式及三级、四级R-K方法实现自动选择步长。L22计算模型为：u=41,（P37,例4.1），精确解为u（t）=（1+t2ju0=1求解：一、PMECME算法1）用四步四阶Adams方法预估，用三步四阶隐式Adams方法校正ph”_一一一一一

2、一i. P：Un+=Un+（55fn-59fn+37fn4-9fn二）24ii. 四步四阶Adams方法的截断误差u（tn+）-Urp+=现卜5/5）tn）+O（h6）（1）720三步四阶隐式Adams方法的截断误差Utn”U：1=-券卜%5tnOh62（1）-（2）得,u：4-Unp+=|0h5uftn。4）两边同乘以-251,则有缉（u：«Unp书）=21h5u（5*tn）+O（h6）720270720iii.iv.v.vi.vii.一一1919两边同乘以，则有9720270(u；甲uc平)=-；170h5u(5Vn)+O(h6)(4)将（3，弋入（1）得,M1：um产u3+25

3、1720E1:fnm1=ftmn：1,un1ch_m._一_一一C：un1-un9fn119fn-5fnJfn_224c19pc将（4y弋入（2）得，M2：un书二un书+商（un中un书E2：fn4=f(%书，un母)编程并作图，将计算解与精确解比较。2）用四步四阶Adams方法预估，用三步四阶Milne方法校正i.ii.hP：un1/2455fn-59忆37%-9人四步四阶Adams方法的截断误差u(tn+)一u：*=H1h5u(5,tn广。化6)三步四阶Milne方法的截断误差uGn+buc+=-Jh5u(5Vn)+O(h6)DU(5)-(6)得,uC-unp+=|9h5uQtn)+O(

4、h6)两边同乘以251259两边同乘以-8-259251则有251259,则有笆h5uiii.iv.v.vi.vii.2597201一905tnOh67h5u(5tn)+O(h6)(8)将（7戕入（5）得，M1：um书二口：4+C：u：1=un-1hfnm14fn+fn43251259uC将（8X弋入（6）得，M2：un由二琮中+:+uc+E2：fn+=f+，4十)编程并作图，将计算解与精确解比较。变步计算1）用误差的后验估计式及三级、四级R-K方法实现自动选择步长。i.三级R-K方法选用Kutta三阶方法（课本P36）,四级R-K方法选用经典的Runge-Kutta法(课本P36,(1.4.

5、22),将两个方法分别编程为函数R-K3.m和R-K4.m。ii.近似解误差的后验估计式unh)-utn：2p2p-1|Lun"iii.iv.要使un(h)-u(tn可利用判断2P又丁>1，最终判断对象为2P-12P5八一0:;2P-12uh-un若“<”成立，则h等于或小于合适的步长，若要进一步判断，可令h=2h,继续下去。若“<”不成立，则h大于合适的步长，将h减半，继续计算直到“<”成立。编写脚本文件zidong.m,不妨令t的取值范围为10,1，固定初始步长为1,对每个t的取值利用ii判断求解合适的步长，并储存所有步长，求出步长的最大和最小值。同时，可

6、储存t的所有取值和对应的计算解，作图与精确解比较以判断误差。程序:PartI编写PMECME算法一clc;clear;h=0.001;%步长t=0:h:1;N=length(t);%区问及划分u(1)=1;f(1)=0;%由Runge-Kutta三阶方法求出初值forj=2:4u(j)=R_K3(t(j),u(j-1),h);f(j)=4*t(j).*sqrt(u(j);endu_p=u;u_c=u;u_m=u;f_m=f;%中间变量赋初值%P%四步四阶%C%三步四阶隐式fori=4:N-1u_p(i+1)=u(i)+h/24*(55*f(i)-59*f(i-1)+37*f(i-2)-9*f(

7、i-3);Adamsu_m(i+1)=u_p(i+1)+251/270*(u_c(i)-u_p(i);%Mf_m(i+1)=4*(i+1)*h*sqrt(u_m(i+1);%Eu_c(i+1)=u(i)+h/24*(9*f_m(i+1)+19*f(i)-5*f(i-1)+f(i-2);Adamsu(i+1)=u_c(i+1)+19/270*(u_p(i+1)-u_c(i+1);%Mf(i+1)=4*(i+1)*h*sqrt(u_m(i+1);%Eendplot(t,u,t,(1+t.A2).A2);legend('PMECME','精确解')PartII编写PM

8、ECME算法二clc;clear;h=0.001;%步长t=0:h:1;N=length(t);%区问及划分u(1:4)=(1+t(1:4)A2)A2;f(1:4)=4*t(1:4).*sqrt(u(1:4);%由精确解求出初值u_p=u;u_c=u;u_m=u;f_m=f;%中间变量赋初值fori=4:N-1u_p(i+1)=u(i)+h/24*(55*f(i)-59*f(i-1)+37*f(i-2)-9*f(i-3);%P%四步四阶Adamsu_m(i+1)=u_p(i+1)+251/259*(u_c(i)-u_p(i);%Mf_m(i+1)=4*(i+1)*h*sqrt(u_m(i+1)

9、;%Eu_c(i+1)=u(i-1)+h/3*(f_m(i+1)+4*f(i)+f(i-1);%C%二步四阶Milneu(i+1)=u_c(i+1)+8/259*(u_p(i+1)-u_c(i+1);%Mf(i+1)=4*(i+1)*h*sqrt(u_m(i+1);%Eendplot(t,u,t,(1+t.A2)A2);legend('PMECME','精确解')Part III R-K3函数functionu=R_K3(t0,u0,h)%Heun三阶方法"k1=4*t0*sqrt(u0);k2=4*(t0+1/3*h)*sqrt(u0+1/3*h*k

10、1);k3=4*(t0+2/3*h)*sqrt(u0+2/3*h*k2);u=u0+h/4*(k1+3*k3);endPart IV R-K4函数functionu=R_K4(t0,u0,h)%四阶Runge-Kutta方法P37(1.4.22)式k1=4*t0*sqrt(u0);k2=4*(t0+1/2*h)*sqrt(u0+1/2*h*k1);k3=4*(t0+1/2*h)*sqrt(u0+1/2*h*k2);k4=4*(t0+h)*sqrt(u0+h*k3);u=u0+h/6*(k1+2*k2+2*k3+k4);endPartV自动变步长clc;clear;H=1;h=H;c=0.000

11、1;%初始步长和精确度t(1)=0;%t的最初取值u(1)=1;%u的估计值k=1;whilet(k)<=1u1=R_K4(t(k),u(k),h);u2=R_K4(t(k),u(k),h/2);whileabs(u1-u2)>=ch=hu1=R_K4(t(k),u(k),h);u2=R_K4(t(k),u(k),h/2);endk=k+1;step(k-1)=h;%合适的步长t(k)=t(k-1)+h;%合适的步长h下t的取值u(k)=u1;%合适的步长h,以h为步长的求解h=H;endext=max(step)min(step);plot(t,u,0:0.001:1,(1+(0:0.001:1)A2)A2);legend('R-K','精确解,);disp('步长的最大值是,num2str(ext(1),'','步长的最小值是,num2str(ext(2)运行结果(1)PMECME算法一(PartI)(2)PMECME算法二(PartII)(3)自动变步长(PartIII、IV、V)(1)