使用C语言解常微分方程CODE.docx

使用C语言解常微分方程CODE.docx解常微分方程
Vincent
年级： 2010，学号： 1033**** ，组号： 5（小组），4（大组）
数值方法 :
我们的实验目标是解常微分方程，其中包括几类问题。一阶常微分初值问题，高阶
常微分初值问题的点。
fk ,1
fk (tm, y1m, y2m,...., ynm )
fk ,2
fk (tm
h , y1m
h f1,1, y2m
h f2,1 ,...., ynm
h f n,1 )
2
2
2
2
fk ,3
fk (tm
h , y1m
h
f1,2 , y2m
h
f2,2 ,...., ynm
h
fn,2 )
2
2
2
2
fk ,2
fk (tm
h, y1m
hf1,3 , y2m
hf2,3 ,...., ynm
hfn,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 的函数
我们对其进行差分
tn
b a
a nh, h
N
是步长
当 h 足够小时，我们有
x ''(tn )
xn 1
2xn
xn 1
h2
x '(tn )
xn 1
xn 1 .
2h
这样的话，我们的微分方程可以写成，
xn 1
2xn
xn 1
h2
p(tn ) xn 1
xn 1
q(tn )xn
r (tn )
0
2h
1
h pn
xn 1
2
h2qn xn
1
h pn xn 1
h2rn , n
1,L ,N 1
2
2
x0 , xN
于是我们得到了个线性方程组
2
h
r
2
h q1
1
2
p1
0
h
2
q2
h
x1
h2 r1
e0
1
2 p2
2 h
1
2 p2
x2
h
2
r2
h pn
h pn
1
2 h2 qn
1
xn
h2 rn
2
2
xN 2
h2rN 2
1
h pN 2
2 h2qN 2
1 h pN 2
xN 1
h2 rN 1
en
r
2
2
1
h p
2 h
2q
0
2
N 1
N
1
其中
e0
( 1 h p1 ) , en
( 1
h pN 1)