使用C语言解常微分方程 C ODE.docVIP

SYSU-IFCEN 2013-2014 实验报告 Projet professionnel Vincent. Wang [PAGE 22] équipe de Neutrons Dosimétrie 解常微分方程 姓名：Vincent 年级：2010，学号：1033****，组号：5（小组），4（大组） 数值方法: 我们的实验目标是解常微分方程，其中包括几类问题。一阶常微分初值问题，高阶常微分初值问题，常微分方程组初值问题，二阶常微分方程边值问题，二阶线性常微分方程边值问题。 对待上面的几类问题，我们分别使用不同的方法。 初值问题 使用 龙格-库塔 来处理 边值问题 用打靶法来处理 线性边值问题 有限差分法 初值问题 我们分别使用 二阶 龙格-库塔 方法 4阶 龙格-库塔 方法 来处理一阶常微分方程。 理论如下： 对于这样一个方程 当h很小时，我们用泰勒展开， 当我们选择正确的参数 a[i][j],b[i][j]之后，就可以近似认为这就是泰勒展开。 对于二阶，我们有： 其中 经过前人的计算机经验，我们有， 选择 A=1/2,B=1/2,则 P=1,Q=1，于是又如下形式，这也使休恩方法的表达式。 所以我们称其为 龙格（库塔）休恩方法 对于4阶龙格方法，我们有类似的想法， 我们使用前人经验的出的系数，有如下公式 对于高阶微分方程及微分方程组 我们用 4阶龙格-库塔方法来解 对于一个如下的微分方程组 我们可以认为是一个一阶向量微分方程，所以可以用龙格-库塔方法的向量形式解。 对于一个高阶的微分方程，形式如下: 我们可以构建出一个一阶的微分方程组， 令 则有 所以我们实际只要解一个微分方程组 其中初值为 使用4阶龙格-库塔方法， 使用这个向量形式的龙格-库塔方法我们便可就出方程的数值解。 边值问题 对于边值问题，我们分为两类 一般的边值问题 线性边值微分方程 一般的边值问题，我们是使用打靶法来求解， 对于这样一个方程 主要思路是，化成初值问题来求解。 我们已有 这样我们便可求出方程的解。 线性微分方程边值问题 对于这样的问题，我们可以使用一个更加高效的方法，有限差分法。 对于如下方程 我们对其进行差分 这样的话，我们的微分方程可以写成， 于是我们得到了个线性方程组 这样的话我们求解出x 对于上面的矩阵，我们可以分解出两个三角阵，然后回代求解便可。 我们用回代法可以直接求解 至此，我们便求出了目标方程的解 流程图 二阶龙格-库塔 对于i = 0到M-1; y[i+1] = y[i] + h / 2 * (fun(t, y[i]) + fun(t + h, y[i] + h*fun(t, y[i]))); return y; 4阶龙格-库塔 对于i = 0到M-1; y[i + 1] = y[i] + h / 6 * (fun(t, y[i]) + 2 * fun(t + h / 2, y[i] + h / 2 * fun(t, y[i])) + 2 * fun(t + h / 2, y[i] + h / 2 * fun(t + h / 2, y[i] + h / 2 * fun(t, y[i]))) +fun(t + h, y[i]+h*fun(t + h / 2, y[i] + h / 2 * fun(t + h / 2, y[i] + h / 2 *fun(t, y[i]))))); return y; 4阶龙格-库塔解方程组 对于k= 0到M-1; 对于i= 0到N; fun(t, y[k], dy[0]) 对于i= 0到N; tempdy[j] = y[k][j] + h / 2 * dy[0][j]; fun(t + h / 2, tempdy, dy[1]); 对于i= 0到N; tempdy[j] = y[k][j] + h / 2 * dy[1][j]; fun(t + h / 2, tempdy, dy[2]); 对于i= 0到N; tempdy[j] = y[k][j] + h * dy[2][j]; fun(t + h, tempdy, dy[3]); y[k+1][i] = y[k][i] + h / 6 * (dy[0][i] + 2 * dy[1][i] + 2 * dy[2][i] + dy[3][i]); return y; 打靶法 当errepsilon y = RKSystem4th( fun, 2, t0, u, tn, M); f0 = y[M-1][0] - b; u[1] = y[0][1]; y = RKSystem4th( fun, 2, t0, v,