常微分方程数值解法上机实验报告.docx

常微分方程数值解法上机实验报告实验题1.ODE:u=1-2tu1+t2u0=0，采用Euler折线法、梯形法、RK3、RK4、三步三阶Adams方法求解。要求：（1）绘制误差下降曲线；（2）确定每个方法的“数值阶”实验图像和数据：（1）误差曲线图：为了比较不同的方法在区间0,T固定、所分区间数N不断变大的情况下，在整个区间上数值解与真解的最大误差下降情况，取T=100、N=1000:10000，得到如下图像：图一：不同方法随着区间数N的增大，最大误差的下降曲线为了探究三步三阶Admas方法的初值选取对后续误差的下降速度的影响，分别用Euler方法给初值，梯形方法给初值，和同阶的RK3方法给初值，同样考虑区间数N变大时，最大误差的下降情况，得到如下图二。 图二：使用不同阶数的方法给出初值对Admas方法的影响（2）数值阶：利用如下方法来计算数值阶：ehT=cheh2T=ch2order=log2eheh2取h=0.1，T=100，对h、h/2、h/4、h/8、h/16，分别计算eh(T)，再利用上述公式，两两相邻计算数值阶，可得如下表格：Euler方法1.0230763111.0111956971.00551595711.002706721梯形方法1.9887729752.0001562161.9984676212.000017355RK33.1010330853.0511392423.0265982583.013430782RK44.0214762184.0112186024.0057134444.002890934三步三阶Admas2.9383190132.98437935832.9960819872.999019688实验数据分析：首先这里函数f(t,u)对u满足 Lipschitz 条件，故初值问题是适定的，且可求得真解为：u(t)=3t+t33+3t2.从图一中我们可以看出，随着区间数的增大，5种方法的误差都在下降，由于问题是适定的，这是合理的。且在图中，高阶方法比低阶的方法下降地更快。我们知道,整体误差e=C*h，其中常数C与函数的高阶导数大小相关，为方法的阶数，故整体误差取决于方法的阶数和高阶导数的大小。在一般情形,即函数的高阶导数增长不剧烈的情况下，则此时高阶方法整体误差较小；但在极端情形下，即函数的高阶导数增长特别剧烈时，也有可能反而较低阶的方法整体误差更小。此处对应函数的高阶导数增长不剧烈，故所得结果合理。从图1中可以看出，一阶方法Euler法的下降最慢，其次三阶Admas方法的误差下降与梯形方法差距不大，再然后RK3和RK4依次下降速度变大。注意到三阶Admas方法的误差下降与梯形方法差距不大，而梯形方法仅为二阶方法，三阶Admas方法与同为3阶的RK3方法相比误差下降速度小了不少。对此做出的解释是：一方面三阶Admas方法的另两个初始值由单步Euler法给出，可能影响到了后序误差的下降速度；另一方面，注意到三步Admas法每多迭代一次只计算一个新的函数值，而RK3要计算3个新的函数值，计算量大了不少，也更多地利用了原问题的信息，从这个角度来看故必然有Admas方法的误差下降速度不如同阶的单步法。.从上面图2中可以看出，对三步三阶Admas 方法，由Euler法、梯形法、RK3给初值的误差下降的效果依次变好，但三者效果几乎平行，并无如图1中Admas方法和RK3所呈现出的那般本质上的区别，且通过不断缩小步长h，三者均能达到同样小的误差界。由此得出结论，不必非得用同阶或更高阶的单步法来给多步法提供初值。.从数值阶表格中可以看出，五种方法的数值阶均与理论上的阶数对应，且注意到虽然在图1中三阶Admas方法的表现效果与RK3方法相比有明显的差距，但两者的数值阶都近似为3，没有大的区别。实验题2.ODE：u=xu(u-2)u0=2，采用如上方法计算；若u0=2+ ，结果如何？实验图像： 设置区间为0,10,初始值为u0=2+10-4,区间数分别为N=100,200,300,500，得到如下数值解与真解的对应图像：同样设置区间为0,10,初始值为u0=2+10-4,步长为0.1,分别使用不同的方法，可得如下图像：实验数据分析：这个初值问题由于对u不满足lipschitz条件，故不是适定的。问题的真解为:u(t)=2u0u0-u0-2et2.取初值u0=2.在此种情形下，问题的真解为ut=2,且此时有ft,2=0, t0.对上述Euler方法、梯形法、RK3、RK4、三步三阶Admas方法对应的迭代公式，从第一步开始每一步都有 fxnyn=0，且利用MATLAB计算时没有受到舍入误差的影响，故最终得到数值解都为yn=2，与真解对应。 取u0=2+. 此时真解为：ut=4+22+-et2.从上面图三可以看出取=10-4时，应用Euler法，当区间数较小，即区间步长较大时，后期数值解会急剧增大。同时可以注意到步长越小，对应的数值解的激增点越靠后。例如：h=110时，从t=3.2开始剧增；而h=130时，从t=9.2才开始剧增。做出的解释如下：对Euler方法的精确计算值为yn+1=yn+h*xnyn(yn-2),由于舍入误差的影响，得到的摄动过的值为，yn+1=yn+h*xnyn(yn-2),令en=yn-yn,则en+1=(1+h*xn(yn+yn- 入误差的影响，得到的摄动过得值555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555552)* en若初值为2，则每一步过后Euler方法的结果仍为2，无误差。若初值为 2+,此时考虑到舍入误差的影响，每一步en被放大(1+h*xn(yn+yn- 入误差的影响，得到的摄动过得值555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555555552)倍，由于初始时刻(yn+yn- 入误差的影响，得到的摄动过得值666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666662)2，当h*xn12 时，放大因子为2，误差开始被急剧放大，则会导致后期误差的剧烈增长。故步长h越小，数值解的激增点越靠后。由于在选择数值方法的步长时，还需考虑到稳定性。步长的变化范围应在绝对稳定区间内。同时注意到Euler法的绝对稳定区间为：(-2,0) ，此处近似地可令=fu=2ux，可以看出任意的步长 h 均不能使 h绝对稳定区间(-2,0),故Euler方法必然难以保证数值解收敛到真解。从上面图四中可以看出，取同一区间0,T，相同步长h=0.1，相同初值u0=2+10-4,五种不同的数值方法给出的数值解的表现也不同。其中由于梯形方法为隐式方法，这里实际上使用的是改进Euler方法。从图中可以观察到：（1）Euler数值解的激增点最靠后；（2）只看激增点之前的时间段，Euler方法的数值解与真解的误差最大，而其他四种方法差距不大，都很好地逼近了真解。上面的五种方法，若取定步长，不管步长h多小，由于相应的=fu=2ux一直在变大 ，故必然无法使h位于绝对稳定区间内。故对于此类问题，必须使用变步长，例如对Euler折线法，基本思路是调整步长h，让每一步的增长因子 (1+h*xn(yn+yn- 入误差的影响，得到的摄动过得值666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666662)不至于太大，若步长小到一定程度，则停机。实验题3.dpdt=p+pqdqdt=q+pq，=-1,=0.01,=0.25,=-0.01初值p0,q0=301,801，绘制相图 (pt,qt）实验图像：下面的图14对应取相同的初值、不同的区间数时，使用Euler折线法时所得的相图。 图1.初值为(31,81) 图2.初值为(29,81) 图3.初值为(31,79) 图4.初值为(29,79) 为了比较不同的初值的影响，下面利用RK4方法，将不同的初值绘制在一张图上，并采用同样的区间数，以便进行比较。图5.区间数为N=200图6.区间数为N=400图7.区间数为N=600实验图像分析：1.从上面图1图4中可以看出，对应同样的初值，使用Euler方法计算时，当步长h较大时，相图趋向于一个闭合的椭圆，当步长h较小时，相图便不再是闭合的椭圆圈了，而是某种不规则图形。2.从上面图5图8中可以看出，对应相同的