机械 CAD程序编制的数学基础编制机械 CAD程序时常用的计算方法
方程求根及其程序
线性方程组求解
数值积分
常微分方程数值解法第一章 方程求根方法方程求根求解区间内的初始近似解牛顿迭代法 弦截法 二分法方程根近似值的确定设方程 f(x)=0,在区间 [a,b]内有一个且只有一个实数根 x*,根据方程的单调性连续性可知，在 x*两侧的函数值 f(x)肯定不同号。根据这个原理，我们便可以确定某个区间 [a,b]内的单实根 x*的近似解。
从左端 x0=a开始，取一个步长 h,h<=b-a,按照步长 h一步一步向区间上限靠近，每增加一个步长进行一次根的扫描，判断 f(x0)与
f(x0+h)是不是异号。如果同号，此时令 x0=x0+h,按步长 h再向右跨。
如果异号，那么 x*必定在区间 [x0,x0+h],可以取 x0或者 x0+h为近似值。
a
b
x *
h
f ( x )
x
注意到 f(0)= -1<0,f(2)=5>0
可见 f(x)在区间 [0,2]之间至少有一个实根。
设从 x=0出发，取 h=0.5为步长向 x=2靠近，记录各个节点上函数的符号，我们发现区间 [1,1.5]
内必有实根，因此可取 x0=1.0或者 x0＝ 1.5作为根的初始近似值。
判断次数 x0 X0+h F(x0) F(x0+h) F(x0)*f(x0+h)
1 0 0.5 -1 -1.375 1.375>0
2 0.5 1 -1.375 -1 1.375>0
3 1 1.5 -1 0.875 -0.875<0
01)( 3 xxxf例题一 确定 的初始近似值。
开始
x0=a
y0=f(x0)
x0=x0+h
输入,a,b,h,f(x)
x0<b f(x0)*y0<=0 输出 x0
初始 x0设定不对 结束
N
N
Y
Y
方程根精确值的确定（ 1）
牛顿迭代法注：牛顿迭代法是先确定根的某个初始近似值，然后用公式反复校正根的近似值，使之逐渐精确化。
)(
)(
'1
k
k
kk xf
xfxx
例题 1：
迭代次数 xk xk+1 xk-xk+1
0 0.5 0.57102 -0.07102
1 0.57102 0.56716 0.00386
2 0.56716 0.56714 -0.00002
弦截法求根,
由于牛顿法需要计算倒数，如果函数 f(x)比较复杂，我们可以使用弦截法，我们可以使用商差来替代牛顿公式中的倒数 f’(x),于是牛顿迭代形式变为：
这个公式的几何意义在于
)()()( )( 1
1
1?
kk
kk
k
kk xxxfxf
xfxx
例题 2：
用弦截法求解方程我们令此时弦截迭代公式为：
01xxe
1)( xxexf
)(
1
)(
)()(
)(
1
1
1
1
1
1
kkx
k
x
k
x
k
k
kk
kk
k
kk
xx
exex
ex
x
xx
xfxf
xf
xx
kk
k
迭代次数 Xk－ 1 xk xk+1
1 0.5 0.6 0.56532
2 0.6 0.56532 0.56715
3 0.56532 0.56715 0.56714
取 x0=0.5,x1=0.6 作为初始近似根二分法求根设已知方程 f(x)=0,且 f(x)在区间 [a,b]内单调连续。显然，若
f(a)*f(b)<0,则方程 f(x)在区间 [a,b]内有且只有一个根 x*。
二分法的基本数学思想是：取方程根所在区间的中点
xm=(a+b)/2，
将所在区间平均分为二个想区间，然后判断根在那个小区间内，
也就是说 f(xm)与 f(a)是否同号。此时有三种情况：
1,f(xm)=0,表示 xm就是方程的根，即 x*=xm;
2,f(xm)f(a)>0,说明根不在 [a,xm]区间内，而是在 [xm,b]区间内。
3,f(xm)f(a)<0,说明根在 [a,xm]区间内。
对有根区间循环上述过程，根的存才范围每次缩小一半，则可以求出一系列根的近似值。在实际计算时候，一般设定一个满足工程需要的计算精度?，使得在有限此次计算之后,|xn-
xn-1|<?，则 xn即所求得近似根。
输入 a,b,?
开始
(a+b)/2=>x
f(a)f(x)<0
If b-a<?
输出 x
结束令 b=x 令 a=x
True
True
false
false
例题三：
使用两分法求方程 在区间（ 1,1.5）内的根。
要求计算精度到 10E-2.
解：由题目知道,a=1,b=1.5,取区间中点
x0=(1+1.5)/2=1.25,由于 f(1.25)= - 0.2969<0与 f(1)同号。则令
a1=1.25,b1=1.5 此时新的有根区间为 [1.25,1.5]。
对新的有根区间取中点 x1=(1.25+1.5)/2=1.375.
由于 f(1.375)=0.2246>0,f(1.25)<0,则令 a2=1.25,b2=1.375，
如此下去就可以得到满足方程精度需求的解。
01)( 3 xxxf
k ak bk xk f(xk)× f(ak)的符号
0 1 1.5 1.25 ＋
1 1.25 1.5 1.375 －
5 1.3125 1.3281 1.3203 ＋
6 1.3203 1.3281 1.3242 ＋
当大区间内有多个根时的求根方法在讨论前面求根方法时，我们是假设方程 f(x)在区间 [a,b]
上 f’(x)及 f”(x)存在且符号不发生改变，这就使得这个函数为单调连续，所以如果 f(a)*f(b)<0,那么就存在唯一的解。但是如果选择区间过大，即使 f(a)*f(b)<0,那么也有可能出现多个解的情况 。
所以，从上面的分析我们可以看出，当区间两端的函数值小于零时，方程至少存在一个根，或者存在奇数个根。当区间两端函数值之积大于零时，方程存在偶数个根或者无根。因此当我们解决大的区间问题时，可以把这个区间分为若干个子区间，如果此子区间足够小且两端函数符号相反，则必定有根，如果符号相同，则无根。
所以当求某函数 f(x)在特定区间的全部根时，过程如下：
1,将区间 [a,b]分成适当大小的 n等分，则区间宽为 h=(b-a)/h
2,从区间下界 a点出发求第一个子区间两端函数值，如果异号则求根，如果同号，则换下个区间，直全部子区间求完。
线性方程组求解,三角形方程组解法基本知识：
线性方程组求解是科学计算中最常遇到的问题，如在应力分析、
电路分析、分子结构、测量学中都会遇到解线性方程组问题,在很多广泛应用的数学问题的数值方法中，如三次样条、最小二乘法、微分方程边值问题的差分法与有限元法也都涉及到求解线性方程组,
n个变量 n个线性方程的方程组我们可以将上述方程组用向量矩阵表示,Ax=b
其中：
其中,A表示 n*n阶实矩阵,x,b为 n维实向量。如果 det(A)不为零，
则方程组有唯一解。理论上来说我们可以用 Cramer法则将方程组求解，但是由于计算量太大，不实用于实际求解。
Cramer法则：
其中?=det(A),而?i表示把?中第 i列的元素换成 b端的对应元素后得到的 n阶行列式。总共需要计算的次数为 N=n!(n*n-1)+n次，对于阶数较大的计算，用此法则无法计算。
nix ii,,2,1,
Gauss法解线性方程组
Gauss消去法就是将方程组通过 (n-1)步消元，转化为上三角方程组，
再回代求此方程组的解,
下面记增广矩阵第一步：设,用第一行加到第 I行，可以消去假定已完成了 (k-1)步消元，即 已将 转化为以下形式：
其中：
照此步骤循环，知道循环 n-1次后，直接回代解
nkkiaala kkkkikkik,,2,1,/,0 1)1()1(11计算 ][ 1)1( bAlik 乘?
),3,2()1( nia ik
][ 1)1( bA
1,.,,,2,1,,.,,,2,1, nkkjnkkiaaaaa
kk
ikkjijij
最初原始方程组变为：
Gauss消去法小节：
1,消元过程对于 k=1,2,…,n-1,若按照顺序有某一 ark不为零,r>=k，
则将 aij-aik*akj/akk=>aij,i=k+1,k+2,…,n; j= k+1,k+2,…,n+1
2,回代过程对于 k=n,n-1,…,,2,1,计算
kk
n
kj
jkjkk axabx /)(
1
三角形分解法（ Doolittle）
原理,对于线性方程组 Ax=b,若其系数矩阵 A能够分解成为两个三角形矩阵相乘,A=LU,其中 L为下三角矩阵,U为上三角矩阵。则方程组 (Ax=b)就变为,LUx=b。 令 Ux=y 则 Ly=b
由图中我们看出：
若已求得 U的 (i-1)行及 L的 (i-1)列，则由矩阵乘法有由上式可得：
所以，关于 L,U的元素可按以下公式计算：
对于 k=1,2,…,,n进行：
1
1
,,1,,
k
s
sjkskjkj nkkjulau?
1
1
,,2,1,/)(
k
s
kkskisikik nkkiuulal?
在确定了两个三角矩阵之后,
数值积分：复合梯形法对于定积分,将积分区间 [a,b]分成 n个相等的子区间，这里步长 h=(b-a)/n＝ xi+1-xi,i=0,1,…,n-1,然后我们在每一个小区间内使用梯形公式：
ba dxxf )(
)]()([2)()( 1
1
0
1
0
1
kkb
a
n
k
n
k
x
x
xfxfhdxxfdxxf k
k
])()(2)([
2
1
1
n
k
k bfxfaf
h
例题：用 n＝ 8的复合梯形公式计算积分解：我们将区间 [0,1]分为 8等份。此时 n=8,h=0.125,f(xk)=sin(xk)/xk
由复合梯形公式可得：
由复合 simpson公式求得的结果为,0.9460832
10 sin dxx x
k xk f(xk) k xk f(xk)
0 0 1.0000000 5 0.625 0.9361556
1 0.125 0.9973978 6 0.75 0.9088516
2 0.25 0.9896158 7 0.875 0.8771926
3 0.375 0.9767267 8 1 0.8414709
4 0.5 0.9588510
9 4 5 6 9 0 9.0])1()(2)0([16 1
7
1
8
k
k fxffT
Simpson法（复合抛物线法）
原理：梯形法用若干段直线组成的折线来替代曲线 f(x),存在误差。如果我们使用若干段抛物线来替代曲线 f(x),误差比梯形法小。
方法：将区间 [a,b]分成
n个相等的子区间
[x2i,x2i+2] (I=0,1,2,3…),设每个子区间上的中点为 x2i+1,且令 h=(b-a)/2n
所以在每一个子区间上利用抛物线公式：
)()(4)(
3
)(
90
)()(4)(
3
)(
22122
)4(
5
22122
12
2
iii
iiii
x
x
xfxfxf
h
f
h
xfxfxf
h
dxxf
i
i
1
1
2
1
12
1
21222
1
21222
)(2)(4)()(
3
)()(4)(
3
)()(4)(
3
n
i
i
n
i
i
n
i
iii
n
i
i
iiii
xfxfbfaf
h
xfxfxf
h
II
xfxfxf
h
I也就是说：将初始区间分成 n个小区间，分别是：
[x0,x2],[x2,x4]等，但是步长 h 是区间的一半。
开 始输 入 a,b,n
( b - a ) / 2 n = > h,a = > x
f ( a ) - f ( b ) = > s
F o r i = 1,n
S = S + s * f ( a + 2 * i * h ) + 4 *
f ( a + 2 * i * h - h )
S = S * h / 3
输 出,返 回开 始输 入 a,b,n
( b - a ) / 2 n = > h,a = > x
f ( a ) - = > s
x + h = > x,
S + 4 f ( x ) = > S
S + f ( b ) = > S
S = S * h / 3
输 出,返 回
b - 2 h < x
x + h = > x
S + 2 f ( x ) = > S
T
F
常微分方程数值解法在工程设计中，常常会遇到常微分方程的求解问题。
这类问题最简单的数学形式试求函数 y(x)在区间 [a,b]上满足一阶常微分方程及其初始条件。
常见的数值解法由很多种，例如 Euler法,Runge-Kutta
法。由于 Euler法计算精度较低，实际应用较多的是
Runge-Kutta法。
00
'
)(
))(,()(
yxy
bxayxfxy
数值解法的基本思想是在区间 [a,b]上取 n+1个离散点，
将初值问题转化成为离散节点的相应问题。然后逐个求出解 yi作为精确解的近似值。
Runge-Kutta法的计算式为：
),(
)
2
,
2
(
)
2
,
2
(
),(
)22(
6
34
23
12
1
43211
hkyhxfk
k
h
y
h
xfk
k
h
y
h
xfk
yxfk
kkkk
h
yy
ii
ii
ii
ii
ii
开 始输 入 a,b,h,y 0
a = > x
y 0 = > y
f ( x,y ) = > k 1
f ( x + h / 2,y + k 1 h / 2 ) = > k 2
f ( x + h / 2,y + k 2 h / 2 ) = > k 3
f ( x + h,y + h k 3 ) = > k 4
y + h ( k 1 + 2 k 2 + 2 k 3 + k 4 ) / 6 = > y
x + h = > x
输 出 x,y
x < b
结 束
T r u e
F a l s e
本章小节
复习和巩固了计算机编程技巧
学习了如何使用计算机进行简单的数学计算，为将来的 CAD设计资料程序化奠定了基础
学习了几种求根，解线性方程组，常微分方程，定积分的常规数值解法上机实践的内容有一台机器，主轴每转一圈，转角?和转矩 T之间的关系如表所示。试计算主轴每转一周所作的功：
1 2 3 4 5 6 7 8 9
Rad
0 0.785 1.571 2.356 3.142 3.927 4.712 5.498 6.283
T
N.mm
5 6.2 7.6 8.1 8.5 9 8.5 6.5 5
510?
求解线性方程组：
68.0
12.1
08.0
76.0
74.008.006.014.0
08.072.004.012.0
06.004.086.002.0
14.012.002.078.0
4
3
2
1
x
x
x
x
方程求根及其程序
线性方程组求解
数值积分
常微分方程数值解法第一章 方程求根方法方程求根求解区间内的初始近似解牛顿迭代法 弦截法 二分法方程根近似值的确定设方程 f(x)=0,在区间 [a,b]内有一个且只有一个实数根 x*,根据方程的单调性连续性可知，在 x*两侧的函数值 f(x)肯定不同号。根据这个原理，我们便可以确定某个区间 [a,b]内的单实根 x*的近似解。
从左端 x0=a开始，取一个步长 h,h<=b-a,按照步长 h一步一步向区间上限靠近，每增加一个步长进行一次根的扫描，判断 f(x0)与
f(x0+h)是不是异号。如果同号，此时令 x0=x0+h,按步长 h再向右跨。
如果异号，那么 x*必定在区间 [x0,x0+h],可以取 x0或者 x0+h为近似值。
a
b
x *
h
f ( x )
x
注意到 f(0)= -1<0,f(2)=5>0
可见 f(x)在区间 [0,2]之间至少有一个实根。
设从 x=0出发，取 h=0.5为步长向 x=2靠近，记录各个节点上函数的符号，我们发现区间 [1,1.5]
内必有实根，因此可取 x0=1.0或者 x0＝ 1.5作为根的初始近似值。
判断次数 x0 X0+h F(x0) F(x0+h) F(x0)*f(x0+h)
1 0 0.5 -1 -1.375 1.375>0
2 0.5 1 -1.375 -1 1.375>0
3 1 1.5 -1 0.875 -0.875<0
01)( 3 xxxf例题一 确定 的初始近似值。
开始
x0=a
y0=f(x0)
x0=x0+h
输入,a,b,h,f(x)
x0<b f(x0)*y0<=0 输出 x0
初始 x0设定不对 结束
N
N
Y
Y
方程根精确值的确定（ 1）
牛顿迭代法注：牛顿迭代法是先确定根的某个初始近似值，然后用公式反复校正根的近似值，使之逐渐精确化。
)(
)(
'1
k
k
kk xf
xfxx
例题 1：
迭代次数 xk xk+1 xk-xk+1
0 0.5 0.57102 -0.07102
1 0.57102 0.56716 0.00386
2 0.56716 0.56714 -0.00002
弦截法求根,
由于牛顿法需要计算倒数，如果函数 f(x)比较复杂，我们可以使用弦截法，我们可以使用商差来替代牛顿公式中的倒数 f’(x),于是牛顿迭代形式变为：
这个公式的几何意义在于
)()()( )( 1
1
1?
kk
kk
k
kk xxxfxf
xfxx
例题 2：
用弦截法求解方程我们令此时弦截迭代公式为：
01xxe
1)( xxexf
)(
1
)(
)()(
)(
1
1
1
1
1
1
kkx
k
x
k
x
k
k
kk
kk
k
kk
xx
exex
ex
x
xx
xfxf
xf
xx
kk
k
迭代次数 Xk－ 1 xk xk+1
1 0.5 0.6 0.56532
2 0.6 0.56532 0.56715
3 0.56532 0.56715 0.56714
取 x0=0.5,x1=0.6 作为初始近似根二分法求根设已知方程 f(x)=0,且 f(x)在区间 [a,b]内单调连续。显然，若
f(a)*f(b)<0,则方程 f(x)在区间 [a,b]内有且只有一个根 x*。
二分法的基本数学思想是：取方程根所在区间的中点
xm=(a+b)/2，
将所在区间平均分为二个想区间，然后判断根在那个小区间内，
也就是说 f(xm)与 f(a)是否同号。此时有三种情况：
1,f(xm)=0,表示 xm就是方程的根，即 x*=xm;
2,f(xm)f(a)>0,说明根不在 [a,xm]区间内，而是在 [xm,b]区间内。
3,f(xm)f(a)<0,说明根在 [a,xm]区间内。
对有根区间循环上述过程，根的存才范围每次缩小一半，则可以求出一系列根的近似值。在实际计算时候，一般设定一个满足工程需要的计算精度?，使得在有限此次计算之后,|xn-
xn-1|<?，则 xn即所求得近似根。
输入 a,b,?
开始
(a+b)/2=>x
f(a)f(x)<0
If b-a<?
输出 x
结束令 b=x 令 a=x
True
True
false
false
例题三：
使用两分法求方程 在区间（ 1,1.5）内的根。
要求计算精度到 10E-2.
解：由题目知道,a=1,b=1.5,取区间中点
x0=(1+1.5)/2=1.25,由于 f(1.25)= - 0.2969<0与 f(1)同号。则令
a1=1.25,b1=1.5 此时新的有根区间为 [1.25,1.5]。
对新的有根区间取中点 x1=(1.25+1.5)/2=1.375.
由于 f(1.375)=0.2246>0,f(1.25)<0,则令 a2=1.25,b2=1.375，
如此下去就可以得到满足方程精度需求的解。
01)( 3 xxxf
k ak bk xk f(xk)× f(ak)的符号
0 1 1.5 1.25 ＋
1 1.25 1.5 1.375 －
5 1.3125 1.3281 1.3203 ＋
6 1.3203 1.3281 1.3242 ＋
当大区间内有多个根时的求根方法在讨论前面求根方法时，我们是假设方程 f(x)在区间 [a,b]
上 f’(x)及 f”(x)存在且符号不发生改变，这就使得这个函数为单调连续，所以如果 f(a)*f(b)<0,那么就存在唯一的解。但是如果选择区间过大，即使 f(a)*f(b)<0,那么也有可能出现多个解的情况 。
所以，从上面的分析我们可以看出，当区间两端的函数值小于零时，方程至少存在一个根，或者存在奇数个根。当区间两端函数值之积大于零时，方程存在偶数个根或者无根。因此当我们解决大的区间问题时，可以把这个区间分为若干个子区间，如果此子区间足够小且两端函数符号相反，则必定有根，如果符号相同，则无根。
所以当求某函数 f(x)在特定区间的全部根时，过程如下：
1,将区间 [a,b]分成适当大小的 n等分，则区间宽为 h=(b-a)/h
2,从区间下界 a点出发求第一个子区间两端函数值，如果异号则求根，如果同号，则换下个区间，直全部子区间求完。
线性方程组求解,三角形方程组解法基本知识：
线性方程组求解是科学计算中最常遇到的问题，如在应力分析、
电路分析、分子结构、测量学中都会遇到解线性方程组问题,在很多广泛应用的数学问题的数值方法中，如三次样条、最小二乘法、微分方程边值问题的差分法与有限元法也都涉及到求解线性方程组,
n个变量 n个线性方程的方程组我们可以将上述方程组用向量矩阵表示,Ax=b
其中：
其中,A表示 n*n阶实矩阵,x,b为 n维实向量。如果 det(A)不为零，
则方程组有唯一解。理论上来说我们可以用 Cramer法则将方程组求解，但是由于计算量太大，不实用于实际求解。
Cramer法则：
其中?=det(A),而?i表示把?中第 i列的元素换成 b端的对应元素后得到的 n阶行列式。总共需要计算的次数为 N=n!(n*n-1)+n次，对于阶数较大的计算，用此法则无法计算。
nix ii,,2,1,
Gauss法解线性方程组
Gauss消去法就是将方程组通过 (n-1)步消元，转化为上三角方程组，
再回代求此方程组的解,
下面记增广矩阵第一步：设,用第一行加到第 I行，可以消去假定已完成了 (k-1)步消元，即 已将 转化为以下形式：
其中：
照此步骤循环，知道循环 n-1次后，直接回代解
nkkiaala kkkkikkik,,2,1,/,0 1)1()1(11计算 ][ 1)1( bAlik 乘?
),3,2()1( nia ik
][ 1)1( bA
1,.,,,2,1,,.,,,2,1, nkkjnkkiaaaaa
kk
ikkjijij
最初原始方程组变为：
Gauss消去法小节：
1,消元过程对于 k=1,2,…,n-1,若按照顺序有某一 ark不为零,r>=k，
则将 aij-aik*akj/akk=>aij,i=k+1,k+2,…,n; j= k+1,k+2,…,n+1
2,回代过程对于 k=n,n-1,…,,2,1,计算
kk
n
kj
jkjkk axabx /)(
1
三角形分解法（ Doolittle）
原理,对于线性方程组 Ax=b,若其系数矩阵 A能够分解成为两个三角形矩阵相乘,A=LU,其中 L为下三角矩阵,U为上三角矩阵。则方程组 (Ax=b)就变为,LUx=b。 令 Ux=y 则 Ly=b
由图中我们看出：
若已求得 U的 (i-1)行及 L的 (i-1)列，则由矩阵乘法有由上式可得：
所以，关于 L,U的元素可按以下公式计算：
对于 k=1,2,…,,n进行：
1
1
,,1,,
k
s
sjkskjkj nkkjulau?
1
1
,,2,1,/)(
k
s
kkskisikik nkkiuulal?
在确定了两个三角矩阵之后,
数值积分：复合梯形法对于定积分,将积分区间 [a,b]分成 n个相等的子区间，这里步长 h=(b-a)/n＝ xi+1-xi,i=0,1,…,n-1,然后我们在每一个小区间内使用梯形公式：
ba dxxf )(
)]()([2)()( 1
1
0
1
0
1
kkb
a
n
k
n
k
x
x
xfxfhdxxfdxxf k
k
])()(2)([
2
1
1
n
k
k bfxfaf
h
例题：用 n＝ 8的复合梯形公式计算积分解：我们将区间 [0,1]分为 8等份。此时 n=8,h=0.125,f(xk)=sin(xk)/xk
由复合梯形公式可得：
由复合 simpson公式求得的结果为,0.9460832
10 sin dxx x
k xk f(xk) k xk f(xk)
0 0 1.0000000 5 0.625 0.9361556
1 0.125 0.9973978 6 0.75 0.9088516
2 0.25 0.9896158 7 0.875 0.8771926
3 0.375 0.9767267 8 1 0.8414709
4 0.5 0.9588510
9 4 5 6 9 0 9.0])1()(2)0([16 1
7
1
8
k
k fxffT
Simpson法（复合抛物线法）
原理：梯形法用若干段直线组成的折线来替代曲线 f(x),存在误差。如果我们使用若干段抛物线来替代曲线 f(x),误差比梯形法小。
方法：将区间 [a,b]分成
n个相等的子区间
[x2i,x2i+2] (I=0,1,2,3…),设每个子区间上的中点为 x2i+1,且令 h=(b-a)/2n
所以在每一个子区间上利用抛物线公式：
)()(4)(
3
)(
90
)()(4)(
3
)(
22122
)4(
5
22122
12
2
iii
iiii
x
x
xfxfxf
h
f
h
xfxfxf
h
dxxf
i
i
1
1
2
1
12
1
21222
1
21222
)(2)(4)()(
3
)()(4)(
3
)()(4)(
3
n
i
i
n
i
i
n
i
iii
n
i
i
iiii
xfxfbfaf
h
xfxfxf
h
II
xfxfxf
h
I也就是说：将初始区间分成 n个小区间，分别是：
[x0,x2],[x2,x4]等，但是步长 h 是区间的一半。
开 始输 入 a,b,n
( b - a ) / 2 n = > h,a = > x
f ( a ) - f ( b ) = > s
F o r i = 1,n
S = S + s * f ( a + 2 * i * h ) + 4 *
f ( a + 2 * i * h - h )
S = S * h / 3
输 出,返 回开 始输 入 a,b,n
( b - a ) / 2 n = > h,a = > x
f ( a ) - = > s
x + h = > x,
S + 4 f ( x ) = > S
S + f ( b ) = > S
S = S * h / 3
输 出,返 回
b - 2 h < x
x + h = > x
S + 2 f ( x ) = > S
T
F
常微分方程数值解法在工程设计中，常常会遇到常微分方程的求解问题。
这类问题最简单的数学形式试求函数 y(x)在区间 [a,b]上满足一阶常微分方程及其初始条件。
常见的数值解法由很多种，例如 Euler法,Runge-Kutta
法。由于 Euler法计算精度较低，实际应用较多的是
Runge-Kutta法。
00
'
)(
))(,()(
yxy
bxayxfxy
数值解法的基本思想是在区间 [a,b]上取 n+1个离散点，
将初值问题转化成为离散节点的相应问题。然后逐个求出解 yi作为精确解的近似值。
Runge-Kutta法的计算式为：
),(
)
2
,
2
(
)
2
,
2
(
),(
)22(
6
34
23
12
1
43211
hkyhxfk
k
h
y
h
xfk
k
h
y
h
xfk
yxfk
kkkk
h
yy
ii
ii
ii
ii
ii
开 始输 入 a,b,h,y 0
a = > x
y 0 = > y
f ( x,y ) = > k 1
f ( x + h / 2,y + k 1 h / 2 ) = > k 2
f ( x + h / 2,y + k 2 h / 2 ) = > k 3
f ( x + h,y + h k 3 ) = > k 4
y + h ( k 1 + 2 k 2 + 2 k 3 + k 4 ) / 6 = > y
x + h = > x
输 出 x,y
x < b
结 束
T r u e
F a l s e
本章小节
复习和巩固了计算机编程技巧
学习了如何使用计算机进行简单的数学计算，为将来的 CAD设计资料程序化奠定了基础
学习了几种求根，解线性方程组，常微分方程，定积分的常规数值解法上机实践的内容有一台机器，主轴每转一圈，转角?和转矩 T之间的关系如表所示。试计算主轴每转一周所作的功：
1 2 3 4 5 6 7 8 9
Rad
0 0.785 1.571 2.356 3.142 3.927 4.712 5.498 6.283
T
N.mm
5 6.2 7.6 8.1 8.5 9 8.5 6.5 5
510?
求解线性方程组：
68.0
12.1
08.0
76.0
74.008.006.014.0
08.072.004.012.0
06.004.086.002.0
14.012.002.078.0
4
3
2
1
x
x
x
x
