xyzsv用c语言编辑,使用C语言解常微分方程CODE.docx

解常微分方程

姓名：Vincent

年级：2010,学号：1033****，组号：5 (小组)，4 (大组)

数值方法：

我们的实验目标是解常微分方程，其中包括几类问题。一阶常微分初值问题，高阶 常微分初值问题，常微分方程组初值问题，二阶常微分方程边值问题，二阶线性常 微分方程边值问题。

对待上面的几类问题，我们分别使用不同的方法。

? 初值问题

使用龙格-库塔来处理

? 边值问题

用打靶法来处理

? 线性边值问题

有限差分法

初值问题

我们分别使用

? 二阶龙格-库塔方法

? 4阶龙格-库塔方法

来处理一阶常微分方程。

理论如下：

对于这样一个方程

y'(t) f(t,y)

当h很小时，我们用泰勒展开，

khf(tk,yQ

k2 hf (tk a〔h,yk KK)

L

i 1

kihf(tk ahM hbjkj)

j 1

当我们选择正确的参数 a[i][j],b[i][j]之后，就可以近似认为这就是泰勒展开。 对于二阶，我们有：

y(t h) y(t) h Af° Bf 1

其中

f。 f (t, y)

f1f(t Ph,y Qhf°)

经过前人的计算机经验，我们有，

选择A=1/2,B=1/2，则P=1,Q=1,于是又如下形式，这也使休恩方法的表达式。

所以我们称其为龙格(库塔)休恩方法

h

y(t h) y(t) f(t,y) f(t h,y hf(t,y)) 2

对于4阶龙格方法，我们有类似的想法， 我们使用前人经验的出的系数，有如下公式

yn 1

yn

h

(k1

6

2k2 L

k1

f(tn,

yn)，

k2

f(tn

h

2 , yn

h

2 k1),

k3

f(tn

h

2,yn

吼), 2

k4

f(tn

h,yn

hk’).

3k4 )>

对于高阶微分方程及微分方程组

我们用

? 4阶龙格-库塔方法来解

对于一个如下的微分方程组 '

y〔f〔(t, y〔,L ,yn),

L

'

ynfn(t,y〔,L , yn).

y〔(t。)ye

L

yn(t°)yn,0

我们可以认为是一个一阶向量微分方程，所以可以用龙格-库塔方法的向量形式解。

TOC \o "1-5" \h \z 对于一个高阶的微分方程，形式如下：

y(n)(x)f t,y,y',L ,y(n1),

\o "Current Document" y(t°)0,

L

\o "Current Document" y(n1)(t°)n1

我们可以构建出一个一阶的微分方程组，

令

y〔(t) y'(t),

L

(n 1)/

yn 〔(t) y (t)

则有

Vn 1(x)f t, y, y〔,L ,yn 1 ,

,

yn 2(x)yn 1(x),

M

,

y2(x)y2(x),

y'(x)y〔(x)

y(to)0,

其中，初值为y1(to)1,

L

V1

f〔(t, y〔,L ,yn),

yq。)

y1,0:

L

其中初值为

L

,

yn

fn(t,y〔,L ,yn).

yn(t°)

yn,0

使用

4阶龙格-库塔方法，

对丁 k §,n ?有,

yn 1(t0)n 1

所以我们实际只要解一个微分方程组

2fk,22fk,3+fk,4

m 1 m h

ykyk 6 fk^

其中,h是步长，

m m

fk (tm, y , y2

m是递归的点。

fk,1

m、

...,yn)

fk,2

fk (tm

h m

2,y1

上f

2 2,1，...

2fn,1)

3, k

? .2

h 一 2myy

,

2

,h 一 2ml

y

,h - 2

m

fk,2fk(tmh,yihfi,3,y2hf2,3,...., ynhfn,3)

使用这个向量形式的龙格-库塔方法我们便可就出方程的数值解。

边值问题

对于边值问题，我们分为两类

? 一般的边值问题

? 线性边值微分方程

一般的边值问题，我们是使用打靶法来求解，

对于这样一个方程

x'' f(t,x, x'), a t b

边界为,x(a) , x(b).

主要思路是，化成初值问题来求解。

我们已有

x(a),

我们估计x'(a) mk

利用初值问题方法求解出x(b) Sk

我们用割线法迭代，使得在进行一定数量的步骤后,

Sk ;

这样我们便可求出方程的解。

?线性微分方程边值问题

对于这样的问题，我们可以使用一个更加高效的方法，有限差分法。 对于如下方程

x''(t) p(t)x'(t) q(t)x(t) r(t) t a,b x(a) , x(b)

其中p,q,r是t的函数

我们对其进行差分

tna nh, h -―-

N

h是步长

当h足够小时，我们有

x''(tn)

x''(tn)

xn 12xn

h^

xn 1

xn 1

2xn