xihy1iK13y2iK23y
iK
3
2程序框图
高阶微分方程转化为一阶微分方程组
取步长
计算K1K2K3K4
yi1
yi
16
K1
K2
K3
K4
3MATLAB编程实现编写4级4阶龙格库塔积分法求解(步长005):
yyyy2x3
y0
1
y0
3
y0
2
转化为一阶常微分方程组:
fy1yy2yy3y
yxy1y2y3T
y0132T
f
x
y
y2
y3
y3
y2
y1
2x
3T
编写MATLAB程序myRKm(见附录),求出总长度为5的数值积分结果,并与解析解作比较。
4算例结果龙格库塔积分结果与实际函数值比较。
8004级RKyxxex2x1
700
600
500
400300
200
100
0
100
0
1
2
3
4
5
6
五、上机总结
通过本次计算方法A课程的上机实验,复习和巩固数值计算方法的基本数学模型,对共轭梯度法求解线性方程组,三次样条插值,牛顿插值,龙贝格积分法,标准四级四阶龙格库塔法求解高阶常微分方程的初值问题等知识有深刻的理解与认识,并通过MATLAB编程实现各种数值算法,对数值分析的应用有清晰的认识。
f附录
1共轭梯度法mygem共轭梯度法fu
ctio
xmygeAble
gthsizeb1x0zerosle
gth1迭代初始向量r0bAx0d0r0fork1le
gth
ar0r0d0Ad0x0x0ad0r1bAx0if
ormr1108
xx0formatlo
gxdisp迭代次数kbreake
dp
ormr12
ormr02d0r1pd0r0r1e
d
2三次样条插值myspli
em三次样条插值
fu
ctio
Sxmyspli
e
a1b1
xiaba
b
yi1125xi2
fori1
hixii1xii
e
d
fori1
1
uihihihi1
vi1ui
e
d
fori2
di16yii1yiihiyiiyii1hi1hihi1
e
d
zuzeros
11zyzeros
11Mzeros
11
M
110M110
第一种边界条件(设为自然边界条件)
d1d1u1M11d
1d
1v
1M
11
fzu112zy11d1追赶法求解三弯矩方程组fori2
1
luizui11zui12lvi1zyi1dilzyi11l0e
dM
1zy
11zu
11forj1
2M
j1zy
j11v
j1M
j11zu
j11e
dMsymsxSxi0set_Sxforj1
Sxixij1x36hjMj1xxij36hjMj11yijhj26Mj1xij1xhjyij1hj26Mj11xxijhj
set_Sxset_SxSxie
dset_SxSxset_Sxforj1
x_Sxxij00001xij1y_Sxsubsset_Sxjx_Sxplotx_Sxy_Sxbholdo
e
d
disp次三次样条插值多项式SxVpaSx4
3牛顿插值myNewto
mmyNewto
mfu
ctio
NxmyNewto
a1b1xiaba
byi1125xi2tablezeros
1
1table1yi差商表tableforj2
1
fforij
1tableijtableij1tablei1j1xiixiij1
e
de
dsymsxNxtable11Tx1fork2
1
TxTxxxik1NxNxr