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

偏微分方程的数值解法 上机习题四 题目： PMECME算法 用四步四阶Adams方法预估（课本P14，k=3），用三步四阶隐式Adams方法校正（课本P16，k=3），编写PMECME算法。 用四步四阶Adams方法预估（课本P14，k=3），用三步四阶Milne方法校正（课本P17，（1.2.31）），编写PMECME算法。 变步计算 用误差的后验估计式及三级、四级R-K方法实现自动选择步长。 计算模型为：，（P37，例4.1 ），精确解为 求解： PMECME算法 用四步四阶Adams方法预估，用三步四阶隐式Adams方法校正 P： ∵四步四阶Adams方法的截断误差 三步四阶隐式Adams方法的截断误差 ∴ 得， 两边同乘以，则有 两边同乘以，则有 将代入得，M1： E1： C： 将代入得，M2： E2： 编程并作图，将计算解与精确解比较。 用四步四阶Adams方法预估，用三步四阶Milne方法校正 P： ∵四步四阶Adams方法的截断误差 三步四阶Milne方法的截断误差 ∴ 得， 两边同乘以，则有 两边同乘以，则有 将代入得，M1： E1： C： 将代入得，M2： E2： 编程并作图，将计算解与精确解比较。 变步计算 用误差的后验估计式及三级、四级R-K方法实现自动选择步长。 三级R-K方法选用Kutta三阶方法（课本P36），四级R-K方法选用经典的Runge-Kutta法（课本P36，（1.4.22）），将两个方法分别编程为函数R-K3.m和R-K4.m。 近似解误差的后验估计式 ， 要使，可利用判断， 又∵，∴最终判断对象为。 若“”成立，则h等于或小于合适的步长，若要进一步判断，可令，继续下去。 若“”不成立，则h大于合适的步长，将h减半，继续计算直到“”成立。 编写脚本文件zidong.m，不妨令t的取值范围为，固定初始步长为1，对每个t的取值利用ii判断求解合适的步长，并储存所有步长，求出步长的最大和最小值。 同时，可储存t的所有取值和对应的计算解，作图与精确解比较以判断误差。 程序： Part I——编写PMECME算法一 clc;clear; h=0.001; %步长 t=0:h:1;N=length(t);% 区间及划分 u(1)=1;f(1)=0;% 由Runge-Kutta三阶方法求出初值 for j=2:4 u(j)=R_K3(t(j),u(j-1),h); f(j)=4*t(j).*sqrt(u(j)); end u_p=u;u_c=u;u_m=u;f_m=f;% 中间变量赋初值 for i=4:N-1 u_p(i+1)=u(i)+h/24*(55*f(i)-59*f(i-1)+37*f(i-2)-9*f(i-3)); %P %四步四阶Adams u_m(i+1)=u_p(i+1)+251/270*(u_c(i)-u_p(i)); %M f_m(i+1)=4*(i+1)*h*sqrt(u_m(i+1)); %E u_c(i+1)=u(i)+h/24*(9*f_m(i+1)+19*f(i)-5*f(i-1)+f(i-2)); %C %三步四阶隐式Adams u(i+1)=u_c(i+1)+19/270*(u_p(i+1)-u_c(i+1)); %M f(i+1)=4*(i+1)*h*sqrt(u_m(i+1)); %E end plot(t,u,t,(1+t.^2).^2); legend(PMECME,精确解) Part II——编写PMECME算法二 clc;clear; h=0.001;%步长 t=0:h:1;N=length(t);% 区间及划分 u(1:4)=(1+t(1:4).^2).^2;f(1:4)=4*t(1:4).*sqrt(u(1:4));%由精确解求出初值u_p=u;u_c=u;u_m=u;f_m=f;% 中间变量赋初值 for i=4:N-1 u_p(i+1)=u(i)+h/24*(55*f(i)-59*f(i-1)+37*f(i-2)-9*f(i-3)); %P %四步四阶Adams u_m(i+1)=u_p(i+1)+251/259*(u_c(i)-u_p(i)); %M f_m(i+1)=4*(i+1)*h*sqrt(u_m(i+1)); %E u_c(i+1)=u(i-1)+h/3*(f_m(i+1)+4*f(i)+f(i-1)); %C %二步四阶Milne u(i+1)=u_c(i+1)+8/259*(u_p(i+1)-u_c(i+1)); %M f(i+1)=4*(i+1)*h*sqrt(u_m(i+1)); %E end plot(t,u,t,(1+t