§5 非线性方程组的解法
设非线性方程组为
n i x x x f n i ,,2,10
),,,(21 == (31)
令 T n T n f f f F x x x x ),,,(,),,,(2121 == 则(31)式记为
F(x)=0 (32)
求x*向量，使得F(x*)≡0的问题，就是非线性方程组求解问题，称x*为方程组（32）式的解向量。
多元非线性方程组与一元非线性方程f(x)=0具有相同的形式，因
此，我们可以类似地讨论它的迭代解法，如不动点迭代法、牛顿迭代法。
一、 简单迭代法（不动点迭代）：
步骤：1.对于非线性方程组: F(x)=0 2.写成等价方程组:x=G(x)
3.构造不动点迭代: x (k+1)
=G(x (k)
)
4.给定初始点x (0)，生成序列{x (k)
} 若
*)(lim
x x k k =∞
→，则x *是方程组的解
程序：
function y=g(x)
y(1)=0.7*sin(x(1))+0.2*cos(x(2)); y(2)=0.7*cos(x(1))-0.2*sin(x(2));
%Iteratepro.m
function y=iteratepro(x) x1=g(x); n=1;
while (norm(x1-x)>=1.0e-6)&(n<=1000) x=x1; x1=g(x); n=n+1; end
x1 n
x0=[0.5 0.5];iteratepro(x0)
二、牛顿迭代法
原理：把解非线性方程的牛顿迭代法，推广到方程组
牛顿迭代的核心是用线性函数近似非线性函数，逐步用线性方程的根近似非线性方程的根。
方法：对方程组F(x)=0，用Taylor 展开方法做线性近似。
设n i x f i ,,2,1),( =具有对),,2,1(n j x j =的二阶偏导数。又设
方程组F(x)=0有近似解),,,()
()(2)(1)(k n k k k x x x x =
将F （x ）在x (k)
点作多元函数的Taylor 展开
)()()()()(k k x x J x F x F -+=+高阶项 (33)
其中矩阵J （Jacobi 矩阵）：为x
（k ）
处i f 的各个一阶偏导数矩阵
????
????
??????????????????????????????=n n n n n n x f x f x f x f x f x
f x f x f x f J
2
1
22212121
11
牛顿迭代过程：
记 )()
(k k x x x
-=? （34）
并略去高阶项，得到近似方程组
)()()(k k x F x J -=? （35）
当det （J ）≠0时，解出)
(k x
?，则有
)1(+k x =)(k x +)(k x ? （36）
再将)
1(+k x
代回（35）式，重复计算，这就是牛顿迭代过程。
牛顿迭代过程的收敛性
??
?=+-==--=0
sin 2.0cos 7.0),(0
cos 2.0sin 7.0),(21y x y y x f y x x y x f
解：第一次迭代： 先求出Jacobi 矩阵
???
? ??+-=y x y x
J cos 2.01sin 7.0sin 2.0cos 7.01 取初值5.0)0()
0(==y x
???
?
?
?=1755.133560
.0095885.038569.0J
det(J)=0.42121≠0,而T y x F )018423.0,011114.0(),()0()
0(-=，由
（35）式，得方程组
???
?
??-=????? ?
?018423.001114.01755.133560.0095885.038569
.0)0(x 解得：,)100145.8,0268234.0(3)
0(T x
-?=?因此得新近似解
50801.0,52682.0)1()
1(==y x
此时
T y x F )1008.2,103.1(),(44)1()1(--??=
第二次迭代：将)
1()
1(,y
x 再代入Jacobi 矩阵中，得方程组
???
?
?????-?-=????? ?
?--44)1(1008.2103.117474.135195.0097288.039491.0x 解得：T x )104689.8,100832.3(54)
1(--?-?-=?，因此迭代两次的近似
解为
50793.0,52651.0)1()1()2()1()1()2(=?+==?+=y y y x x x
T y x F )104,0(),(6)2()2(-?=，精度已相当高。
牛顿迭代计算机算法
步1，输入初始解向量ε,)
0(X
；
步2，对i,j=1,2,…,n 计算
)
0(x x f i
i ??； 步3，求)()
0(x
F ；
步4，求解方程组);()0()
0(x F x J -=?
步5，求)0()0()
1(x x x ?+=；
步6，若
ε<-)0()1(max i i i
x x ，则输出)1(x ，转步7否则
)1()0(x x ?，转步2；
步7，结束。 附牛顿迭代法程序： function y=fun5(x)
y(1)=x(1)-0.7*sin(x(1))-0.2*cos(x(2)); y(2)=x(2)-0.7*cos(x(1))+0.2*sin(x(2)); y=[y(1) y(2)];
function y=df(x)
y=[1-0.7*cos(x(1)) 0.2*sin(x(2))
0. 7*sin(x(1)) 1+0.2*cos(x(2))];
%newtonpro.m
function y=newtonpro(x0) x1=x0-fun5(x0)/df(x0); n=1;
while(norm(x1-x0)>=1.0e-6)&(n<=1000) x0=x1;
x1=x0-fun5(x0)/df(x0);n=n+1; end
x1,n %输出满足精度的迭代次数与解 x0=[0.5 0.5];newtonpro(x0)
三、最速下降法
问题的提出：对于方程组 F(x)=0 ，构造各方程的平方和函数，即令
∑==Φn
i i n x f x x x 1
221)(),,,( （37）
则求F(x)=0 的解等价于求多元函数),,,(21n x x x Φ的极小值问题。 这种无约束极值问题，可以采用最优化方法求解。
最速下降法(梯度法):1847年由Cauchy 提出，它是求解最优化问题的最古老且基本的数值方法。后来提出的一些好的方法、都是这种方法的改进。
原理：设),,,(21n x x x Φ对各变元具有连续的偏导数，记其梯度为
T
n
x x x x ),,,(
)(21?Φ??Φ??Φ?=Φ? （38） 取T
n x x x x
),,,()0()0(2)0(1)
0( =为近似解，任取一个单位向量p 及正数t ，
将)(x Φ在)
0(x
点进行Taylor 展开，即有
)()(()()()0()0()0(t o p x t x tp x T +Φ?+Φ=+Φ （39）
其中o(t)是t 的高阶无穷小量，而
)
),(cos()()
),(cos()()]([)
0()0()0()0()0(p x x p x p x p x T Φ??Φ?=Φ???Φ?=Φ? （40）
1.选方向：当（p x
),()
0(Φ?）=π时，即
,)(/)()0()0(时x x p Φ?Φ-?=)()0(tp x +Φ向)()0(x Φ变化最快，就是
说沿负梯度方向)(x Φ在点)
0(x 附近变化最快。最速下降法正是因此而得
名。
2．选步长：从)
0(x
求)
1(x
，沿负梯度方向跨出一步：t ，即
)()0()0()1(x t x x Φ?-= （41）
选步长原则：使))(()()0()
0(x t x
t F Φ?-Φ=达极小，即求0t 使
)(min ))(()0(0)0(t F x t x =Φ?-Φ （42）
重复上述过程，由)
1(x
求出新的近似解)()1(1)1()
2(x t x x
Φ?-=
最速下降法算法公式：k=0,1,2,…
?????Φ?-=Φ?-Φ=Φ?-Φ+)
44()()43())((min ))(()()()1()
()()()(k k k k k k i k k k x t x x
x t x x t x
迭代计算，直至Φ的值降到很小，则终止计算。
【说明】最速下降法对任选的初始向量都能收敛，但收敛速度不如牛
顿迭代法快。
将牛顿迭代法与最速下降法联合使用，可取长补短，效果更好。 最优化原理的使用：详见Matlab 优化工具箱
Matlab 中非线性方程求根问题的实现：
fsolve :Solve a system of nonlinear equations by a least squares method.
【例】function F = fun7(x)
F = [2*x(1) - x(2) - exp(-x(1)); -x(1) + 2*x(2) - exp(-x(2))];
x0=[-5; -5]; % Make a starting guess at the solution options=optimset('Display','iter');
% Option to display output
[x,fval] = fsolve(@fun7,x0,options) % Call optimizer %x=fsolve(@fun7,x0)
x =
0.5671
0.5671
【例】fun = inline('sin(3*x)');
x = fsolve(fun,[1
4],optimset('Display','off'));x
x =
1.0472 4.1888
fun = inline('sin(3*x)');
x = fsolve(fun,1,optimset('Display','off'));x
x =
1.0472
【例】解方程组：x-0.7*sin(x)-0.2*cos(y)=0
y-0.7*cos(x)+0.2*sin(y)=0
x0=[0 0];
fsolve('fun5',x0)
ans =
0.5265 0.5079
solve:symbolic solution of algebraic equations
【例】solve('sin(x) = 1')
ans =
1/2*pi
【例】 solve('p*sin(x) = r') %chooses 'x' as the unknown
ans =
asin(r/p)
【例】方程组[x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0')
x =
[ 1]
[ 3]
y =
[ 1]
[ -3/2]
【例】
[x,y]=solve('x-0.7*sin(x)-0.2*cos(y)','y-
0.7*cos(x)+0.2*sin(y)')
x =
.52652262191818418730769280519209
y =
.50791971903684924497183722688768
第八章 矩阵特征值问题计算 3．用幂法计算下列矩阵的主特征值及对应的特征向量 12732343()341;()463213331a A b A --???? ????=-=-???? ????--???? 当特征值有3位小数稳定代终止。 解：套用幂法公式 010,,,1,2,.... max()k k k k k v u v Au u k v -≠== = 取0(1,1,1)0T u =≠,将A 1代入上式，计算结果见下表 则1A 的主特征值19.605572λ≈，特征向量1(10.6050.394369)T x ≈- 将2A 代入幂法公式，取0(1,1,1)T u =，计算结果见下表 则2A 主住特征值18.869699λ≈，特征向量1(0.604228,1,0.160881)T x ≈- 4．用反幂法求矩阵 621231111A ?? ??=?? ???? 的最接近于6的特征向量。 解：本题按带原点平移的反幂法计算。平移向量p=6,则将
021231115B A pI ?? ??=-=-?? ??-?? 进行三分解：PB=LU ，其中 1 002310101511 001,10,02 221004 2701005 5P L U ? ??? ????-??? ??? ??????===-???????????? ?? ?? ??? ??? 然后1(1,1,1)T Uv =，解得 1 111,max()v v u v = 1,,,2,3,.... max()k k k k k k k v Ly PU Uv y U k v -=== = 计算结果如下：
非线性方程数值解法及其应用 摘要:数值计算方法主要研究如何运用计算机去获得数学问题的数值解的理论和算法。 本文主要介绍非线性方程的数值解法以及它在各个领域的应用。是直接从方程出发，逐步缩小根的存在区间，或逐步将根的近似值精确化，直到满足问题对精度的要求。我将从二分法、Steffensen 加速收敛法、Newton 迭代法、弦截法来分析非线性方程的解法及应用。 关键字：非线性方程;二分法；Steffensen 加速收敛法；代数Newton 法；弦截法 一、前言 随着科技技术的飞速发展，科学计算越来越显示出其重要性。科学计算的应用之广已遍及各行各业，例如气象资料的分析图像，飞机、汽车及轮船的外形设计，高科技研究等都离不开科学计算。因此经常需要求非线性方程 f(x) = O 的根。方程f(x) = O 的根叫做函数f(x)的零点。由连续函数的特性知：若f(x)在闭区间[a ，b]上连续，且f(a)·f(b) 数值分析第8章 数值积分与数值微分 8.1 填空题 （1）n+1个点的插值型数值积分公式∫f(x)dx b a ≈∑A j n j=0f(x j )的代数精度至少是 n ，最高不超过 2n+1 。【注：第1空，见定理8.1】 （2）梯形公式有 1 次代数精度，Simpson 公司有 3 次代数精度。【注：分别见定理8.1,8.3】 （3）求积公式∫f(x)dx h 0≈h 2[f (0)+f (h )]+ah 2[f ′(0)?f ′(h)]中的参数a= 1/12 时，才能保证该求积公式的代数精度达到最高，最高代数精度为 3 。 解：令f(x)=1,x,x 2带入有， { h 2[1+1]+ah 2[0?0]=h h 2[0+h ]+ah 2[1?1]=12 (h 2)h 2[0+h 2]+ah 2[0?2h ]=13 (h 3) //注：x 的导数=1 解之得，a=1/12，此时求积公式至少具有2次代数精度。 ∴ 积分公式为：∫f(x)dx h 0≈h 2[f (0)+f (h )]+h 2 12[f ′(0)?f ′(h)] 令 f(x)= x 3带入求积公式有：h 2 [0 +h 3]+ h 212 [0?3h 2]=14 (h 4)，与f(x)= x 4的定积分计算值1 4 (h 4)相等， 所以，此求积公式至少具有3次代数精度。 令f(x)= x 4带入求积公式有，h 2[0+h 4]+h 2 12[0?4h 3]=1 6(h 5)，与f(x)= x 5的定积分计算值1 5(h 5)不相等，所以，此求积公式的最高代数精度为3次代数精度。 8.2 确定下列求积公式的求积系数和求积节点，使其代数精度尽量高，并指出其最高代数精度。 解题思路：按照P149 中8.3式进行求解，根据求积公式中未知量n 的数量决定代入多少f(x)，当积分公式代入求积节点x n 的计算结果与定积分的计算结果一致，继续代入求积节点X n+1,，若计算结果与对应的定积分计算结果不一致时，求积公式拥有最高n 次的代数精度。 （1）∫f(x)dx 2h 0≈A 0f (0)+A 1f (h )+A 2f(2h) 解：令f(x)=1,x,x 2代入有,【注：本例中需求解A 0、A 1、A 2共3个未知量，故需3个相异求积节点f(x)】 {A 0+A 1+A 2=2h A 1h +A 22h =1 2(2h )2A 1h 2+A 2(2h )2=1 3(2h )3 求解得A 0=13h ，A 1=43h ，A 2=1 3h ， ∴求积公式为：∫f(x)dx 2h 0≈13hf (0)+43hf (h )+1 3 hf(2h) ∵该求积公式对3个相异节点1,x,x 2均有余项E (f )=0， //注：参见P149定理8.1 ∴该求积公式至少具有2次代数精度。 令f(x)= x 3，代入求积公式有：4 3hh 3+1 3h (2h )3=4h 4 ∵函数f(x) = x 3的定积分结果为：∫x 3dx 2h 0=1 4(2h )4=4h 4 ，与求积公式计算值相等， ∴该求积公式具有3次代数精度。 数值分析讲义 第三章线性方程组的解法 §3.0 引言 §3.1 雅可比(Jacobi)迭代法 §3.2 高斯-塞德尔(Gauss-Seidel)迭代法 §3.3 超松驰迭代法§3.7 三角分解法 §3.4 迭代法的收敛性§3.8 追赶法 §3.5 高斯消去法§3.9 其它应用 §3.6 高斯主元素消去法§3.10 误差分析 §3 作业讲评3 §3.11 总结 §3.0 引言 重要性：解线性代数方程组的有效方法在计算数学和科学计算中具有特殊的地位和作用.如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题. 分类：线性方程组的解法可分为直接法和迭代法两种方法. (a) 直接法：对于给定的方程组，在没有舍入误差的假设下，能在预定的运算次数内求得精确解.最基本的直接法是Gauss消去法，重要的直接法全都受到Gauss消去法的启发.计算代价高. (b) 迭代法：基于一定的递推格式，产生逼近方程组精确解的近似序列.收敛性是其为迭代法的前提，此外，存在收敛速度与误差估计问题.简单实用,诱人. §3.1 雅可比Jacobi 迭代法 (AX =b ) 1 基本思想: 与解f (x )=0 的不动点迭代相类似,将AX =b 改写为X =BX +f 的形式,建立雅可比方法的迭代格式：X k +1=BX (k )+f ,其中，B 称为迭代矩阵.其计算精度可控,特别适用于求解系数为大型稀疏矩阵(sparse matrices)的方程组. 2 问题: (a) 如何建立迭代格式？ (b) 向量序列{X k }是否收敛以及收敛条件? 3 例题分析: 考虑解方程组??? ??=+--=-+-=--2.453.82102 .72103 21321321x x x x x x x x x (1) 其准确解为X *={1, 1.2, 1.3}. 建立与式(1)相等价的形式： ??? ??++=++=++=84.02.01.083.02.01.072 .02.01.02 13312321x x x x x x x x x (2) 据此建立迭代公式： ?????++=++=++=+++84 .02.01.083.02.01.072.02.01.0)(2)(1)1(3 )(3 )(1)1(23)(2)1(1k k k k k k k k k x x x x x x x x x (3) 取迭代初值0) 0(3 )0(2)0(1===x x x ,迭代结果如下表. JocabiMethodP31.cpp 第六章非线性方程的数值解法习题解答 填空题： 1. 求方程()x f x =根的牛顿迭代格式是__________________。 Ans:1()1()n n n n n x f x x x f x +-=- '- 2.求解方程 在(1, 2)内根的下列迭代法中， (1) (2) (3) (4) 收敛的迭代法是（A ）. A ．(1)和(2) B. (2)和(3) C. (3)和(4) D. (4)和(1) 3.若0)()(，故迭代发散。 以上三中以第二种迭代格式较好。 2、设方程()0f x =有根，且'0()m f x M <≤≤。试证明由迭代格式1()k k k x x f x λ+=- (0,1,2,)k =L 产生的迭代序列{}0k k x ∞ =对任意的初值0(,)x ∈-∞+∞，当2 0M λ<< 时，均收敛于方程的根。 实验报告 一、实验目的 1．迭代函数对收敛性的影响。 2．初值的选择对收敛性的影响。 二、实验题目 1．用简单迭代法求方程01)(3=--=x x x f 的根。 分别化方程为如下等价方程： 31+=x x ；13 -=x x ；x x 11+=；213-+=x x x 取初值5.10=x ，精度为4 10-，最大迭代次数为500，观察其计算结果并加以分析。 2．①用牛顿法求方程01)(3=-+=x x x f 在0.5附近的根， 分别取初值1000,100,2,1,5.0,5.0,1,2,100,10000-----=x 观察并比较计算结果，并加以分析。 ②用牛顿法求方程0)(3=-=x x x f 所有根。 三、实验原理 简单迭代法程序，牛顿迭代法程序。 四、实验内容及结果 五、实验结果分析 (1)实验1中用简单迭代法求方程01)(3=--=x x x f 的根： 取初始值5.10=x 的时候，等价方程2和4是不收敛的。等价方程1的迭代次数为6，近似值为1.324719474534364。等价方程3的迭代次数为7，近似值为1.324718688942791。说明不同的等价方程得到的结果以及迭代的次数是不一样的。 (2)实验2中用牛顿迭代法求方程01)(3=-+=x x x f 在0.5附近的根： 通过结果可知，当初始值越接近真实值时，迭代的次数就越少。 (3)实验3中用牛顿法求方程0)(3=-=x x x f 所有根： 可知该方程的根为01=x ，12=x ，13-=x ，由于方程是无重根的，所以可以直接用牛顿迭代法做，而不需要使用牛顿迭代加速法做。 一 实验目的 1．掌握求解线性方程组的高斯消元法及列主元素法； 2. 掌握求解线性方程组的克劳特法； 3. 掌握求解线性方程组的平方根法。 二 实验内容 1．用高斯消元法求解方程组（精度要求为610-=ε）： 1231231 233272212240x x x x x x x x x -+=??-+-=-??-+=? 2．用克劳特法求解上述方程组（精度要求为610-=ε）。 3. 用平方根法求解上述方程组（精度要求为610-=ε）。 4. 用列主元素法求解方程组（精度要求为610-=ε）： 1231231 233432222325x x x x x x x x x -+=??-+-=??--=-? 三 实验步骤（算法）与结果 1. 程序代码（Python3.6）： import numpy as np def Gauss(A,b): n=len(b) for i in range(n-1): if A[i,i]!=0: for j in range(i+1,n): m=-A[j,i]/A[i,i] A[j,i:n]=A[j,i:n]+m*A[i,i:n] b[j]=b[j]+m*b[i] for k in range(n-1,-1,-1): b[k]=(b[k]-sum(A[k,(k+1):n]*b[(k+1):n]))/A[k,k] print(b) 运行函数： >>> A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=np.float) >>> b=np.array([7,-1,0],dtype=np.float) >>> x=Gauss(A,b) 输出： 结果：解得原方程的解为x1=3.5,x2=-1,x3=-2.25 2 程序代码（Python3.6）： import numpy as np A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=float) L=np.array([[1,0,0],[0,1,0],[0,0,1]],dtype=float) U=np.array([[0,0,0],[0,0,0],[0,0,0]],dtype=float) b=np.array([7,-1,0],dtype=float) y=np.array([0,0,0],dtype=float) x=np.array([0,0,0],dtype=float) def LU(A): n=len(A[0]) i=0 while i scanf("%lf",&X0); printf("请输入精度：\n"); scanf("%lf",&tol); printf("请输入最大循环次数\n"); scanf("%d",&n); X1=X0-f(X0)/ft(X0); while(fabs(X1-X0)>=tol)/*在没有达到精度要求情况下就循环*/ { X0=X1; X1=X0-f(X0)/ft(X0); count++;/*循环次数计数*/ /*如果循环次数到用户输入的最大循环次数就退出，可以适当修改*/ if(count==n) { break; } } printf("方程的根为:%f\n",X1); getch(); } /************************* 函数功能：用于求f(x)的值 **************************/ double f(double x) { double y; y=x*x*x+2*x*x+10*x-20; return y; } /************************* 函数功能：用于求f(x)的导数值 **************************/ doubleft(double x) { double y; y=3*x*x+4*x+10; return y; } 四、例题 用Newton 迭代法求下面方程的一个正根,计算结果精确到7位小数. 02010223=-++x x x []得取初值,2,020∈=x 郑州大学研究生课程（2012-2013学年第一学期）数值分析 Numerical Analysis 习题课 第八章常微分方程数值解法 待求解的问题：一阶常微分方程的初值问题/* Initial-Value Problem */: ?????=∈=0 )(] ,[),(y a y b a x y x f dx dy 解的存在唯一性（“常微分方程”理论）：只要f (x , y ) 在[a , b ] ×R 1 上连续，且关于y 满足Lipschitz 条件，即存在与x , y 无关的常数L 使 对任意定义在[a , b ] 上的y 1(x ) 和y 2(x ) 都成立，则上述IVP 存在唯一解。 1212|(,)(,)||| f x y f x y L y y ?≤?一、要点回顾 §8.2 欧拉(Euler)法 通常取(常数),则Euler 法的计算格式 h h x x i i i ==?+1?? ?=+=+) (),(001x y y y x hf y y i i i i i =0,1,…,n ( 8.2 ) §8.2 欧拉(Euler)法(1) 用差商近似导数 )) (,()()()()(1n n n n n n x y x hf x y x y h x y x y +=′+≈+?? ?=+=+) (),(01a y y y x hf y y n n n n 差分方程初值问题向前Euler 方法h x y x y x y n n n ) ()()(1?≈ ′+)) (,() ()(1n n n n x y x f h x y x y ≈?+))(,()(n n n x y x f x y =′ 第八章习题解答 3、设方程()0f x =有根，且'0()m f x M <≤≤。试证明由迭代格式1()k k k x x f x λ+=- (0,1,2,)k =产生的迭代序列{}0k k x ∞=对任意的初值0(,)x ∈-∞+∞，当20M λ<<时，均收敛于方程的根。 证明： 设()()x x f x ?λ=-，可知()x ?在(,)-∞∞上可导 对于任意给定的λ值，满足条件'0()m f x M <≤≤时 （1）''()1()x f x ?λ=- 则1'()11M x m λ?λ-≤≤-< 又20M λ<<，M>0 则02M λ<<时，11M λ-<- 所以11'()11M x m λ?λ-<-≤≤-< 若令max{1,1}L M m λλ=--，则可知'()1x L ?≤< （2）由0()(0)'()(0)'()x x x dx x ?????ε=+=+? 则()lim 1x x L x ?→∞??≤< ??? 所以，存在一个数a ，当x a >时，()x x ?< 同时，()x ?在[,]a a -内有界，即存在0b >使得[,]x a a ?∈-，()x b ?< 我们选取 max{,}c a b =，则对任意x 有0()max{,}x c x ?< 则对给定的任意初值0x ，设0max{,}d c x = 则0[,]x d d ∈-，于是在区间[,]d d -上有()x d ?< 即满足映内性 有（1）、（2）可知，()x ?满足收敛定理 迭代序列0{}k k x ∞=收敛于方程的根 6. 给出计算...222+++=x 的迭代格式，讨论迭代格式的收敛性，并证明2=x 解：构造迭代格式10,1,2,k x k +==??? 2k x ≤ 令()x ?=x ?∈?时，()x ??∈? '() x ?=，当x ?∈?时，1 '()12x ?<< 文章编号:1671-8755(2003)03-0050-05 收稿日期:2002-03-03 四川省教育厅基金项目(2002A022) 非线性数值分析在岩土工程中的应用 苏华友 (西南科技大学环境工程学院 四川绵阳 621010) 摘要:天然岩土体中都存在着不同成因的节理、裂隙和断层等软弱面,在对这种天然材料进行理论 研究与工程设计时,需要应用非线性数值分析方法。介绍了几种主要的非线性数值理论,并重点 讨论了有限元非线性分析和位移反分析法在边坡工程和地下工程中的应用。 关键词:非线性数值分析 位移反分析 岩土工程 中图分类号:O29;TB11文献标识码:A 随着我国经济的飞速发展,许多基础性的产业得到了大力开发和加强,特别是我国西部大开发方针更为我们提供了新的挑战和机遇,交通、水利水电、矿产资源以及许多公共工程显然处于较为突出的地位,这些领域的开发所依附的母体正是我们脚下的岩土。 早期人们对岩土的变形性状、破坏机制以及所建立的力学模型都仅限于线弹性范围内。但是,随着人们对岩石作为一种天然材料与其赋存的地质环境、形成历史,地壳运动以及工程因素之间联系的认识逐步加深,特别是一些对大型土木工程破坏的深入研究,如:法国的马尔巴赛(M alpasset)双曲拱坝,左坝肩边坡部分岩体1959年12月2日突然溃决,香港中心区一座27楼塌滑到山脚等,人们认识到把岩土体简单地当作线弹性材料来处理是远远不够的,需要对岩土体进行非线性数值分析。目前应用较多的非线性数值方法有:弹塑性有限元数值分析法、边界元法、离散元法、断裂力学、位移反分析、损伤力学分析等。 1 几种主要的非线性数值分析法简介 1.1 弹塑性理论求解岩石力学问题 在岩石力学中使用弹塑性理论是将岩石介质看做是一种连续介质(当然它们可以含有有限多个间断面),因为应力、应变等概念都是建立在连续介质模型基础上的。岩石是一种杂质体,它由各种不同的矿物组成,存在节理裂隙等不连续面和空隙,在细观上看,也是不连续的。弹塑性理论抓住问题的主要方面,在宏观的小尺度范围内来考虑各种力学量的统计平均值。 弹塑性力学的研究内容可分为以下两大部分:(1)研究材料固有的特性,建立应力、应变及温度之间关系的本构关系(数字表达式)。(2)分析弹塑性变形物体内应力、应变分布,研究物体在各种载荷作用下的稳定性问题。 1.2 边界元法 边界元法是同有限元法并行发展的另一类数值法。该方法在岩石力学中的应用自20世纪70年代以后有了较大的发展。P.K.Banerjee,C.A.Brebbia,S.L.Crouch 等人都曾对其发展和应用作出过有价值的贡献。我国自20世纪70年代进行岩石力学边界元的应用研究,在岩体稳定性分析、地下工程支护等方面做出了有第18卷 第3期 西南科技大学学报 Vol.18No.3 2003年9月JOU RN AL O F SOU TH WEST U N IV ERSIT Y OF SCIENCE A ND T ECHN OLO GY Sep.2003 非线性回归分析 1．问题描述 在熔盐泵模化试验中，根据模化方案，采用水作为试验介质，对输送介质密度ρp =1938kg/m3，粘度μp=0.00729Pa?s，转速n p=1450r/min，颗粒直径1mm，密度1938kg/m3，叶轮直径D p=250mm的原型泵进行模化试验，取x=0，得模型泵的转速n m=386.2r/min，叶轮直径D m=250mm。颗粒1mm，直径密度1103.8kg/m3。 在不同的流量工况下，分别对这两台泵进行数值模拟，然后将流量、扬程和轴功率全部转换成无量纲的比流量、比扬程和比功率。对其中的效率，比能量数值曲线进行拟合，实验数据结果如下，其中效率为a，比能量为b。 试求效率y对比能量x的回归方程。 方法一： 2. 用STATISTICA进行非线性回归分析 根据y与x的对应数据，EXCEL绘图可以看出来，他们之间满足指数关系（如下图所示），所以设回归方程为b-Bb2+Cb3。用STATISTICA做试验分析时采用自定义回归方程模块。 建立如图所示的数据表： 3.回归过程详解 输入自定义回归方程： 采用Levenberg-Marquardt估计方法求解结果显示对话框： 分析结果显示： 方差分析结果如下： 观测值，预测值，残差值如下： 残差直方图： 残差散点图： 观测值与回归曲线对比图： 结果分析： 从上面的分析结果里我们可以看到系数a=18.0639，b=-91.9099，c=147.391，d=-0.42801即y=18.0639x3-91.9099x2+147.391x-0.42801。我们可以看出拟合曲线和散点之间的相关度是0.99984089。从残差直方图可以看出图像近似满足正态分布规律。从残差散点图可以看出残差点没有明显的规律可寻，即说明残差基本满足随机分布。综合以上分析可以说明，预测的回归曲线方程的参数与解析解非常的接近。关系曲线方程为y=18.0639x3-91.9099x2+147.391x-0.42801。 方法二： 使用origin进行非线性回归： 1.原曲线： 数值分析习题集 （适合课程《数值方法A》和《数值方法B》） 长沙理工大学 第一章绪论 1.设x>0,x的相对误差为δ,求的误差. 2.设x的相对误差为2％,求的相对误差. 3.下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指 出它们是几位有效数字: 4.利用公式求下列各近似值的误差限: 其中均为第3题所给的数. 5.计算球体积要使相对误差限为1％,问度量半径R时允许的相对误差限是多少? 6.设按递推公式 ( n=1,2,…) 计算到.若取≈(五位有效数字),试问计算将有多大误差? 7.求方程的两个根,使它至少具有四位有效数字(≈. 8.当N充分大时,怎样求? 9.正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝? 10.设假定g是准确的,而对t的测量有±秒的误差,证明当t增加时S的绝对误差增加,而 相对误差却减小. 11.序列满足递推关系(n=1,2,…),若(三位有效数字),计算到时误差有多大?这个计算过程 稳定吗? 12.计算,取,利用下列等式计算,哪一个得到的结果最好? 13.,求f(30)的值.若开平方用六位函数表,问求对数时误差有多大?若改用另一等价公式 计算,求对数时误差有多大? 14.试用消元法解方程组假定只用三位数计算,问结果是否可靠? 15.已知三角形面积其中c为弧度,,且测量a ,b ,c的误差分别为证明面积的误差满足 第二章插值法 1.根据定义的范德蒙行列式,令 证明是n次多项式,它的根是,且 . 2.当x= 1 , -1 , 2 时, f(x)= 0 , -3 , 4 ,求f(x)的二次插值多项式. 3. 4.给出cos x,0°≤x ≤90°的函数表,步长h =1′=(1/60)°,若函数表具有5位有效数 字,研究用线性插值求cos x 近似值时的总误差界. 5.设,k=0,1,2,3,求. 6.设为互异节点(j=0,1,…,n),求证: i) ii) 7.设且,求证 8.在上给出的等距节点函数表,若用二次插值求的近似值,要使截断误差不超过,问使用函 数表的步长应取多少? 9.若,求及. 10.如果是次多项式,记,证明的阶差分是次多项式,并且为正整数). 11.证明. 12.证明 13.证明 14.若有个不同实根,证明 15.证明阶均差有下列性质: i)若,则; ii)若,则. 16.,求及. 17.证明两点三次埃尔米特插值余项是 并由此求出分段三次埃尔米特插值的误差限. 18.求一个次数不高于4次的多项式,使它满足并由此求出分段三次埃尔米特插值的误差限. 19.试求出一个最高次数不高于4次的函数多项式,以便使它能够满足以下边界条件,,. 20.设,把分为等分,试构造一个台阶形的零次分段插值函数并证明当时,在上一致收敛到. 21.设,在上取,按等距节点求分段线性插值函数,计算各节点间中点处的与的值,并估计误 差. 22.求在上的分段线性插值函数,并估计误差. 23.求在上的分段埃尔米特插值,并估计误差. i) ii) 25.若,是三次样条函数,证明 i); ii)若,式中为插值节点,且,则. 26.编出计算三次样条函数系数及其在插值节点中点的值的程序框图(可用式的表达式). 第三章函数逼近与计算 1.(a)利用区间变换推出区间为的伯恩斯坦多项式. (b)对在上求1次和三次伯恩斯坦多项式并画出图形,并与相应的马克劳林级数部分和误 实验报告一：实验题目 一、 实验目的 掌握求解非线性方程根的二分法、简单迭代法和牛顿迭代法，并通过数值实验比较两种方法的收敛速度。 二、 实验内容 1、编写二分法、牛顿迭代法程序，并使用这两个程序计算 02)(=-+=x e x x f 在[0, 1]区间的解，要求误差小于 4 10- ，比较两种方法收敛速度。 2、在利率问题中，若贷款额为20万元，月还款额为2160元，还期为10年，则年利率为多少？请使用牛顿迭代法求解。 3、由中子迁移理论，燃料棒的临界长度为下面方程的根cot x =(x 2?1)/2x ，用牛顿迭代法求这个方程的最小正根。 4、用牛顿法求方程f (x )=x 3?11x 2+32x ?28=0的根，精确至8位有效数字。比较牛顿迭代法算单根和重根的收敛速度，并用改进的牛顿迭代法计算重根。 三、 实验程序 第1题： 02)(=-+=x e x x f 区间[0,1] 函数画图可得函数零点约为0.5。 画图函数： function Test1() % f(x) 示意图, f(x) = x + exp(x) - 2; f(x) = 0 r = 0:0.01:1; y = r + exp(r) - 2 plot(r, y); grid on 二分法程序： 计算调用函数：[c,num]=bisect(0,1,1e-4) function [c,num]=bisect(a,b,delta) %Input –a,b 是取值区间范围 % -delta 是允许误差 %Output -c 牛顿迭代法最后计算所得零点值 % -num 是迭代次数 ya = a + exp(a) - 2; yb = b + exp(b) - 2; if ya * yb>0 return; end for k=1:100 c=(a+b)/2; yc= c + exp(c) - 2; if abs(yc)<=delta a=c; b=c; elseif yb*yc>0 b=c; yb=yc; else a=c; ya=yc; end if abs(b-a) 第八章 最优化问题 最优化分支：线性规划，整数规划，几何规划，非线性规划，动态规划。又称规划论。 应用最优化方法解决问题时一般有以下几个特点： 1. 实用性强 2. 采用定量分析的科学手段 3. 计算量大，必须借助于计算机 4. 理论涉及面广 应用领域：工业，农业，交通运输，能源开发，经济计划，企业 管理，军事作战……。 §8.1 最优化问题实例 最优化问题：追求最优目标的数学问题。 经典最优化理论： （1） 无约束极值问题：),,,(opt 21n x x x f （),,,(min 21n x x x f 或),,,(max 21n x x x f ） 其中，),,,(21n x x x f 是定义在n 维空间上的可微函数。 解法（求极值点）：求驻点，即满足 ????? ? ?='='='0 ),,(0),,(0 ),,(1 1121n x n x n x x x f x x f x x f n 并验证这些驻点是否极值点。 （2） 约束极值问题：),,,(opt 21n x x x f s.t. )(,,2,1,0),,,(21n l l j x x x h n j <== 解法：采用Lagrange 乘子法，即将问题转化为求Lagrange 函数 ) ,,(),,,(),,;,,,(11 21121n j j l j n l n x x h x x x f x x x L λλλ∑=+=的无约束极值问题。 近代最优化理论的实例： 例1 (生产计划问题) 设某工厂有3种资源B 1，B 2，B 3，数量各为b 1，b 2，b 3，要生产10种产品A 1，…，A 10 。每生产一个单位的A j 需要消耗B i 的量为a ij ，根据合同规定，产品A j 的量不少于d j ，再设A j 的单价为c j 。问如何安排生产计划，才能既完成合同，又使总收入最多？（线性规划问题） 数学模型：设A j 的计划产量为 j x ，z 为总收入。 目标函数： max 10 1 j j j x c z ∑== 约束条件： ?????=≥=≤∑=10 ,,2,1,3 ,2,1,10 1 j d x i b x a j j j i j ij 线性规划问题通常采用单纯形法来求解。 例2 (工厂设址问题) 要在m 个不同地点计划修建m 个规模不完全相同的工厂，他们的生产能力分别是m a a a ,,21 （为简便起见，假设生产同一种产品），第i 个工厂的建设费用m i f i ,,2,1, =。又有n 个零售商店销售这种产品，对这种产品的需求量分别为 数值分析中求解非线性方程的MATLAB求解程序（6种） 1.求解不动点 function [k,p,err,P]=fixpt(g,p0,tol,max1) %求解方程x=g(x) 的近似值，初始值为p0 %迭代式为Pn+1=g(Pn) %迭代条件为：在迭代范围内满足|k|<1(根及附近且包含初值)k为斜率 P(1)=p0; for k=2:max1 P(k)=feval(g,P(k-1)); err=abs(P(k)-P(k-1)); relerr=err/(abs(P(k))+eps); p=P(k); if (err if b-a 实验名称： 实验五 线性方程组的数值解法 指导教师： 数值分析实验组 实验时数： 2 实验设备：安装了Matlab 、C ++、VF 软件的计算机 实验日期：2014年 月 日 实验地点： 第五教学楼北802或902 实验目的： 1. 掌握线性方程组的迭代法和直接法的基本思想和基本步骤。 2. 理解各类数值解法的优缺点，并能自行编程求解。 3. 认识迭代法收敛的含义以及迭代法初值和方程组系数矩阵性质对收敛速度的影响，了解求解病态线性方程组的方法。 实验准备： 1. 在开始本实验之前，请回顾教科书的相关内容； 2. 需要一台准备安装Windows XP Professional 操作系统和装有数学软件的计算机。 实验内容及要求 A 题 考虑方程组b Hx =的求解，其中系数矩阵H 为Hilbert 阵， n j i j i h h H j i n n j i ,,2,1,,1 1 ,)(,, =-+= =? 这是一个著名的病态问题。通过首先给定解（例如取各个分量均为1）再计算出右端的办法给出确定的问题。 （1） 选择问题的维数为6，分别用列主元Gauss 消去法、Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？ （2） 逐步增大问题的维数，仍然用上述方法来解它们，计算的结果如何？计算的结果说明了什么？分析产生结果的原因。 说明：实验过程应包括对问题的简要分析、求解方法、求解步骤、程序及其必要的图表等内容。 实验过程： 本实验所选题为A 题 实验分析：b Hx =，H 矩阵可由Matlab 直接给出，为了设定参考解，首先设x 为分量全 为1的向量，求出b ，然后将H 和b 作为已知量，求x ，与设定的参考解对比。 对于列主元Gauss 消去法，Jacobi 迭代，Gauss-Seidel 迭代，SOR 迭代法，去迭取迭代初值 第六章非线性方程的数值解法习题解答 填空题： 1. 求方程()x f x =根的牛顿迭代格式是__________________。 Ans:1()1()n n n n n x f x x x f x +-=- '- 2.求解方程 在(1, 2)内根的下列迭代法中， (1) (2) (3) (4) 收敛的迭代法是（A ）. A ．(1)和(2) B. (2)和(3) C. (3)和(4) D. (4)和(1) 3.若0)()(，故迭代发散。 以上三中以第二种迭代格式较好。 2、设方程()0f x =有根，且'0()m f x M <≤≤。试证明由迭代格式1()k k k x x f x λ+=- ）《数值分析》杨大地-标准答案(第八章)
数值分析讲义线性方程组的解法
第六章非线性方程的数值解法习题解答
数值分析非线性方程求根实验
数值分析实验2_求解线性方程组直接法
郑州大学研究生课程数值分析复习---第八章 常微分方程数值解法
最新(完美版)第八章习题答案_数值分析
非线性数值分析在岩土工程中的应用
数值分析-非线性回归分析
数值分析习题集及答案
数值分析求解非线性方程根的二分法,简单迭代法和牛顿迭代法
数值计算方法第八章
数值分析中求解非线性方程的MATLAB求解程序(6种)
数值分析— 线性方程组的数值解法
第六章非线性方程的数值解法习题解答
相关主题
文本预览