数值分析常微分方程数值解实验题.docVIP

数值分析第四次上机练习实验报告 ——常微分方程数值解实验题 问题的描述 考虑一下问题： ； 初试条件为，。其精确解为 ； 请分别用古典四级四阶显式Runge-Kutta方法和隐式二级四阶Runge-Kutta方法计算，计算区间取成，并与精确解比较。 方法描述 Runge-Kutta方法是采用不同点上函数值的不同组合来提高方法的精度。又避免了函数f的偏导数计算。 古典四级四阶显式Runge-Kutta方法 古典四级四阶显式Runge-Kutta方法相应的计算格式为 隐式二级四阶Runge-Kutta方法 隐式二级四阶Runge-Kutta方法相应的计算格式为： 方案设计 我们通过编写程序来进行运算，程序语言采用Matlab语言，运行环境为Matlab R2010b。 在我们的程序文件中practice4.m文件为主程序文件， RK_ex44.m 为使用古典四级四阶显式Runge-Kutta方法计算的函数，RK_im24.m 为使用隐式二级四阶Runge-Kutta方法计算的函数。 古典四级四阶显式Runge-Kutta方法 由function y=RK_ex44(f,a,b,y0,N,p)函数计算得到。该函数可计算函数f，计算区间上下限为a、b，初值为y0，N个点，p个未知数，使用古典四级四阶显式Runge-Kutta方法计算得到的值。 显式计算方法非常便于实现，按照计算公式使用matlab语句编程实现即可。 具体实现详见源程序。 隐式二级四阶Runge-Kutta方法 由function y=RK_im24(A,c,a,b,y0,N) 函数计算得到。该函数可计算函数f，计算区间上下限为a、b，初值为y0，N个点使用隐式二级四阶Runge-Kutta方法计算得到的值。 隐式计算方法不似显式那么便利，需要解非线性方程组。然后依然是迭代求解。 具体实现详见源程序。 计算结果及误差分析 如上图所示为步长h=0.001时的仿真结果。图中蓝色线代表精确解，红色线代表古典显式四级四阶Runge-Kutta方法的解，绿色线（由于颜色叠加显示为黄色线）代表隐式二级四阶Runge-Kutta方法的解。由图可知，当步长h选择足够小时，两种方法均能较好的对实验结果进行模拟。 下面对当选择不同的步长h时的实验结果进行讨论。首先讨论古典显式四级四阶Runge-Kutta方法。 下述六幅图分别显示了当步长h=0.0015，h=0.0014，h=0.0013时，u(t)和v(t)的仿真结果。由图中可以看出，当步长不够小时，无法得到正确的仿真结果。而一旦选择足够小的步长，能够得到较好的仿真结果。  当h=0.0015时 当h=0.0014时 当h=0.0013时 当h=0.1时 接下来讨论步长h对隐式二级四阶Runge-Kutta方法的仿真结果的影响。 当h=0.01时 可以看出，当h=0.1时，隐式二级四阶算法已能较好的精确v(t)，然而，因为u(t)值出现了阶跃，故并不能很好的对u(t)进行仿真。 下表显示了h=0.1,h=0.08, h=0.06, h=0.04, h=0.02, 0.01时，u(t)的仿真图。并根据局部讨论的需要进行了局部的放大。由图中可以看出，随着步长的逐步缩减，在t=0处由于，u(t)的阶跃造成的误差在逐步缩减，直至可以忽略。 结论 通过实验可以发现，显式RK方法虽然表达方式简明，步骤简单，然而需要较小的步长以及较长的计算时间才能达到较精确的解。相比之下，隐式RK方法需要的计算时间远小于显式方法，并且较大的步长就可以满足得到精确解的需求，亦即隐式方法的绝对稳定性较好。不足之处在于，隐式方法表达公式不够简明，需要求解非线性方程组。 在实际生活中，可以根据所需求解问题的需要选取相应的求解方法。