常微分方程初值问题的稳定性研究与实例分析

1实验2常微分方程初值问题的稳定性研究与实例分析一、实验目的主要研究在广泛应用领域中进行科学工程计算时经常出现的常微分方程数值求解的稳定性问题，理解设计数值求解算法不仅要求具有收敛性还需特别重视稳定性。二、实验内容1、从理论上讨论常微分方程解析解关于初值的稳定性及其与常微分方程数值解关于初值的稳定性的不同和关系、总结线性扰动分析理论分析常微分方程数值解关于初值的稳定性的关键步骤；2、应用线性扰动分析理论对自选的几个常微分方程数值求解算法进行稳定性分析讨论3、选择示例进行算法设计，用Matlab进行程序设计与运行实验。三、实验原理、方法和手段3.1线性扰动分析理论对于m阶齐次常系数差分方程∑𝑎𝑘𝑦𝑛+𝑘𝑚𝑘=0=0其中𝑎𝑘为常数，利用移位算子E可定义特征方程为𝑓(𝜆)=∑𝑎𝑘𝜆𝑘𝑚𝑘=0=0可证明差分方程稳定的充要条件是|𝜆𝑖|1对差分格式，令𝐴=𝑓𝑦，则对y给个小扰动δy,则可解得δy=𝑒𝐴(𝑥−𝑥0)𝛿𝑦02当A0时误差随x的增加指数上升，故我们总在A0的条件下讨论稳定性。3.2单步法求解常微分方程初值问题的单步法是指以给定的初值y(x0)和步长h为基础，一步一步计算y(x1),y(x2),…,y(xn)。常用的单步法有欧拉(Euler)方法和龙格-库塔(Runge-Kutta)法，其导出都是依据泰勒（Taylor）展开。迭代格式如下。Euler法：𝑦𝑛+1=𝑦𝑛+ℎ𝑓(𝑥𝑛,𝑦𝑛)Runge-Kutta法：𝑦𝑛+1=𝑦𝑛∑𝜇𝑗𝐾𝑗𝑚𝑗=1𝐾𝑗=ℎ𝑓(𝑥𝑛+𝑎𝑗ℎ,𝑦𝑛+∑𝑏𝑗𝑙𝐾𝑙𝑗−1𝑙=1)𝑗=1,2,…,𝑚3.3多步法单步法中的Euler法计算量小但精度低，Runge-Kutta法正相反。提出多步法的目的是寻求精度高，运算量少的计算公式。多步法的原理是已求得前m个解的前提下，可使用已知的所有数值解来求第m+1个解，这样就必然出现积分问题。在区间[𝑥𝑛,𝑥𝑛+1]对初值问题积分得𝑦(𝑥𝑛+1)=𝑦(𝑥𝑛)+∫𝑓(𝑥,𝑦(𝑥))𝑑𝑥𝑥𝑛+1𝑥𝑛对右边的积分使用拉格朗日差值并舍去余项得m+1歩开型阿达姆斯(Adams)公式𝑦𝑛+1=𝑦𝑛+ℎ∑𝑎𝑘(𝑚)𝑓(𝑥𝑛−𝑘,𝑦𝑛−𝑘)𝑚𝑘=0其中3𝑎𝑘(𝑚)=(−1)𝑘𝑘!(𝑚−𝑘)!∫∏(𝑡+𝑗)𝑑𝑡𝑚𝑗=0𝑗≠𝑘10,𝑘=0,1,…,𝑚四、实验条件台式计算机，Matlab软件试验平台。五、实验步骤5.1线性扰动分析根据第三节介绍的线性扰动分析理论对单步法的Euler格式，Runge-Kutta格式以及多步法的m=3的开型阿达姆斯公式进行扰动分析。(1)Euler格式𝑦𝑛+1=𝑦𝑛+ℎ𝑓(𝑥𝑛,𝑦𝑛)将𝑓(𝑥𝑛,𝑦𝑛)=𝐴𝑦𝑛代入得特征方程𝜑(𝜆)=𝜆−(1+ℎ𝐴)=0得𝜆=(1+ℎ𝐴)。由于A0,h0,故𝜆1。当hA-2时Euler格式是稳定的。即ℎ2|𝐴|(2)Runge-Kutta格式{𝐾1=ℎ𝑓(𝑥𝑛,𝑦𝑛)𝐾2=ℎ𝑓(𝑥𝑛+12ℎ,𝑦𝑛+12𝐾1)𝐾3=ℎ𝑓(𝑥𝑛+12ℎ,𝑦𝑛+12𝐾2)𝐾4=ℎ𝑓(𝑥𝑛+ℎ,𝑦𝑛+𝐾3)𝑦𝑛+1=𝑦𝑛+16(𝐾1+2𝐾2+2𝐾3+𝐾4)此时用线性固化方法将𝑓(𝑥𝑛,𝑦𝑛)=𝐴𝑦𝑛代入得λ=124((ℎ𝐴)4+4(ℎ𝐴)3+12(ℎ𝐴)2+24ℎ𝐴+24)由于hA无限趋于负无穷时λ≈124(ℎ𝐴)4必然能大于1，而Ha=-1时λ0,故只有4hA满足一定条件使|λ|1才能使其稳定。实际大概的稳定条件为ℎ2.7|𝐴|(3)m=3的开型Adams格式𝑦𝑛+1=𝑦𝑛+124ℎ[55𝑓(𝑥𝑛,𝑦𝑛)−59𝑓(𝑥𝑛−1,𝑦𝑛−1)+37𝑓(𝑥𝑛−2,𝑦𝑛−2)−9𝑓(𝑥𝑛−3,𝑦𝑛−3)]将𝑓(𝑥𝑛,𝑦𝑛)=𝐴𝑦𝑛代入得特征方程𝜆4−(1+5524ℎ𝐴)𝜆3+5924ℎ𝐴𝜆2−3724ℎ𝐴𝜆+924ℎ𝐴=0分析可得稳定的条件为ℎ0.3|𝐴|5.2算法设计Euler算法：①初始化，h=(b-a)/N;x=a;y=y0;②对n=1,2,…,N循环y=y+hf(x,y);x=x+h;③输出x,y.Runge-Kutta算法：①初始化h=(b-a)/N;x=a;y=y0;②对n=1,2,…,N循环K1=hf(x,y);K2=hf(x+0.5h,y+0.5K1);5K3=hf(x+0.5h,y+0.5K2);K4=hf(x+h,y+K3);y=y+(K1+2K2+2K3+K4)/6;x=x+h;③输出x,y。m=3的开型Adams算法①初始化，h=(b-a)/N;x=a;𝑦0=𝑦0;②用单步法计算𝑦1,𝑦2,𝑦3.③对n=3,4,…,N-1循环𝑦𝑛+1=𝑦𝑛+124ℎ[55𝑓(𝑥𝑛,𝑦𝑛)−59𝑓(𝑥𝑛−1,𝑦𝑛−1)+37𝑓(𝑥𝑛−2,𝑦𝑛−2)−9𝑓(𝑥𝑛−3,𝑦𝑛−3)]𝑥𝑛+1=𝑥𝑛+ℎ④输出x,y5.3实例选取与编程实现选取课本P102页例4.8.𝑑𝑦𝑑𝑥=−𝑦+𝑥+1,0≤𝑥≤30,𝑦(0)=1解析解为y=𝑒−𝑥+𝑥,此例中A=-1，由5.1的分析可知三种格式的稳定条件分别为Euler:h2;Runge-Kutta:h2.7;Adams:h0.3.考虑到迭代次数必须为整数，可分别取h=0.2,0.3,0.4,1,1.88,2,2.14,2.5,2.73,3。6对b=30处的误差e做如下处理𝑒={0,𝑖𝑓(𝑒10−6)∞,𝑖𝑓(𝑒10)𝑒,𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒若0e∞,且迭代过程中误差基本保持不变，则在数字e后加C作为标记。则实验结果如表5-1所示。其中eE,eR,eA分别表示Euler法，Runge-Kutta法和Adams法在终点b处的误差。表5-1三种迭代格式在终点b处的误差e随步长变化规律步长h0.20.30.41.01.882.02.142.52.733.0次数N1501007530161514121110误差eE00000.1181.0C6.48∞∞∞误差eR0000002×10−60.00550.38C∞误差eA09.2×10−5𝐶∞∞∞∞∞∞∞∞5.4结果分析由表5-1可得，Euler法在h2.0时是稳定的；h=2.0时迭代过程中误差基本保持不变，可认为处于临界稳定状态；h2.0时，迭代过程中误差增大，不具稳定性。同样可得Runge-Kutta法的稳定条件为h2.73,Adams法的稳定条件为h0.3.这样的结果是和5.1节的理论分析相符的。六、实验总结通过本实验，从理论上分析了常系数微分方程初值问题的几种迭代格式的稳定性，实例研究结果与理论分析相符。7附件2实验2代码%filename:DifferentialSolution.m%常微分初值问题%选取例子为课本的P102例4.8%df/dx=-y+x+10=x=30y(0)=1clear;clc;a=0;b=30;y0=1;h=2;N=round(b/h);%单步法的欧拉算法Euler(a,b,y0,N);%单步法的龙格-库塔算法RungeKutta(a,b,y0,N);%多步法的m=3的开型阿达姆斯算法Adams(a,b,y0,N);%filename:Euler.m%欧拉Euler法functionEuler(a,b,y0,N)x=[];y=[];y(1,1)=y0;h=(b-a)/N;fori=1:Nx(i)=a+(i-1)*h;y(i+1,1)=y(i)+h*f(x(i),y(i,1));y(i,2)=fr(x(i));endy(N+1,2)=fr(b);y(:,3)=abs(y(:,2)-y(:,1));%输出disp('Euler方法');fprintf('步长h=%f\n',h);disp('yy_reale');%disp(y);fprintf('%2.10f%2.10f%2.10f\n',y');%filename:RungeKutta.m%四阶Runge-Kutta方法8functionRungeKutta(a,b,y0,N)x=[];y=[];x(1)=a;y(1,1)=y0;h=(b-a)/N;fori=1:NK1=h*f(x(i),y(i,1));K2=h*f(x(i)+h/2,y(i,1)+K1/2);K3=h*f(x(i)+h/2,y(i,1)+K2/2);K4=h*f(x(i)+h,y(i,1)+K3);y(i+1,1)=y(i,1)+(K1+2*K2+2*K3+K4)/6;y(i,2)=fr(x(i));x(i+1)=a+i*h;endy(N+1,2)=fr(b);y(:,3)=abs(y(:,2)-y(:,1));%格式输出disp('四阶Runge-Kutta法');fprintf('步长h=%f\n',h);disp('yy_reale');fprintf('%2.10f%2.10f%2.10f\n',y');%filename:Adams.m%m=3的开型Adams法functionAdams(a,b,y0,N)x=[];y=[];x(1)=a;y(1,1)=y0;h=(b-a)/N;fori=1:3x(i+1)=a+i*h;y(i+1,1)=fr(x(i+1));y(i,2)=fr(x(i));endfori=4:Nx(i)=a+(i-1)*h;y(i+1,1)=y(i)+h*(55*f(x(i),y(i,1))+(-59)*f(x(i-1),y(i-1,1))+37*f(x(i-2),y(i-2,1))+(-9)*f(x(i-3),y(i-3,1)))/24;y(i,2)=fr(x(i));endy(N+1,2)=fr(b);9y(:,3)=abs(y(:,2)-y(:,1));%格式输出disp('开型Adams法');fprintf('步长h=%f\n',h);disp('yy_reale');fprintf('%2.10f%2.10f%2.10f\n',y');%filename:f.m%实例的格式function[fv]=f(x,y)fv=-y+x+1;end%filename:fr.m%实例的解析解function[f]=fr(x)f=exp(-x)+x;end