微分方程c语言求解,使用C语言解常微分方程 C ODE.doc

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;

打靶法

当err

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,