非线性方程求解

disp([’k=’,num2str(k)])
牛顿法
牛顿法也称为牛顿-拉普森法或者切线法。由 于这个方法的计算结果颇佳，而计算过程也比 较简单，所以被普遍采用。
牛顿法的核心内容是通过泰勒级数将非线性方 程式转化为线性方程式，然后用迭代法求解。
牛顿法原理
设方程式 f x 0 的近似根为 x0
数值计算的误差
模型误差
误 差 观测误差
的
分 类
截断误差
计算时只截 ln(x 1) x 1 x2 1 x3 1 x4 (1)n1 1 xn
取有限项
234
ቤተ መጻሕፍቲ ባይዱ
n
舍入误差
计算机对所 储存的数据 位数有限制
误差的分类 有效数字
绝对误差 相对误差
(x) x*x
(x)
二分法 不动点迭代 牛顿法
二分法
若函数f(x)在区间[a,b]内单调连续，且f(a)f(b)<0， 则在闭区间[a,b]内必然存在方程f(x)=0的根x*
二分法的图形解释
二分法的MATLAB程序
function x=mabisec(fun,a,b,ep)
x=(a+b)/2.0; k=0;
课堂练习
某蒸馏釜的操作压力位106.7kpa，其中溶液含苯摩尔分 数为0.2，甲苯为0.8，求此溶液的泡点和平衡的气相组 成。
苯甲苯溶液可以看做理想溶液，组分蒸汽压为：a为苯， b为甲苯
log Pa 6.301 1211 t 220.8
log Pb 6.080 1345 t 219.5
{ x0, x1, x2, …, xk, … }
如果这个数列有极限，则迭代格式是收敛的。
这时数列{xk}的极限
x*

lim
k
xk
就是方程的根
上述求非线性代数方程式数值解的方法称为直
接迭代法（或称为不动点迭代法）。这个方法
虽然简单，但根本问题在于当k->∞时，xk是否
收敛于x*，也就是必须找出收敛的充分条件
得迭代格式如下：
1
xk1 2
xk2 3 ，k
0,1,2,
求得x1

1 2
x02 3
6.5
发散的迭代 格式
x2
1 2
x12 3
19.625
x3 191.070,
当k越来越大时，xk的值越来越远离于精确根
不动点
定义：函数g(x)的一个不动点（fixed point） 是指一个实数x*，满足x*= g(x*)
避免两个相近的数相减 避免用绝对值对过小的数作除数 避免两个数量级相差很大的数运算 算法中应尽量减少运算次数 递推运算中，防止误差的积累
非线性方程求解
非线性方程数值求解基本原理 Matlab求解非线性方程函数－roots， fzero， fsolve
引言
在945.36kPa（9.33atm）、300.2K时，容器中 充以2mol氮气，试求容器体积。已知此状态下 氮气的P-V-T关系符合范德华方程，其范德华常 数为a＝4.17atm•L/mol2， b＝0.0371L/mol
从图形角度分析，函数y=g(x)的不动点是 y=g(x)和y=x的交点
不动点定理
设函数P(x)在区间[a,b]上有连续的一阶导数，且满 足：
a P(x) b, x(a,b); p'(x) L1, x(a,b)
则函数P(x)在区间[a,b]上存在唯一的不动点X*，即有 x*=P(x*)，且对任意的属于区间[a,b] 的x0，都可以通 过迭代得到X*。
while k<N x=x1-(x1-x0)*feval(fun,x1)/(feval(fun,x1)-feval(fun,x0));
if abs(x-x1)<ep break; end x0=x1; x1=x; k=k+1; end if k==N, warning(’ 已达迭代次数上限’); end disp([’k=’,num2str(k)])
正确的收敛 的迭代格式
x3 3.034, x4 3.011, x5 3.004, 当k越来越大时，xk的值越来越接近于精确根
x* 3
例题：f x x2 2x 3 0
如果构造迭代格式如下：
1）x 1 2
x2 3 ,即x 1
2
x2 3 ，取x0 4,
对于多项式，用以下行向量表示：
P [a0 , a1,, an1, an ]
这样就把多项式问题转化为向量问题 Matlab提供了多种多项式计算函数，如多项式求根函数
roots，求多项式的值，polyval；多项式乘法，conv；多 项式除法，deconv；多项式微分，polyder；多项式拟合， polyfit

f xk f ' xk
牛顿法的几何意义
Y
y f x
x0, f x0 x1, f x1
切线方程
y f (x0 ) (x x0 ) f '(x0 )
X
O
x*
x1
x0
例:牛顿法计算x^2-25=0的解
f(x)=x2-25，则f’(x)=2x
可构造迭代公式如下：
牛顿迭代法收敛速度快，但它要求计算函数导数的值
MATLAB求解非线性方程方法
MATLAB求解非线性方程函数
非线性方程
多项式函数
非线性方程
非线性方程组
roots
fzero
fsolve
多项式求根函数roots
多项式的表达式约定如下：
P(x) ao x n a1 x n1 an1 x an
设x是数x的近似值如果x的绝对误差限是他的某一位数位的半个单位并且从x左起第一个非零数字到该数位共有n位则称这个n个数字为x的有效数字也称x近似x时有n位有效数字有效数字算法设计的原则避免两个数量级相差很大的数运算算法中应尽量减少运算次数避免两个相近的数相减递推运算中防止误差的积累避免用绝对值对过小的数作除数matlab求解非线性方程函数rootsfzerofsolve非线性方程数值求解基本原理引言引言在94536kpa933atm3002k时容器中充以2mol氮气试求容器体积
r(x) (x*x) / x r(x) (x*x) / x*
r(x)
定义：设X*是数X的近似值，如果X*的绝对误差限时他 的某一位数位的半个单位，并且从X*左起第一个非零数 字到该数位共有N为，则称这个N个数字为X*的有效数字， 也称X*近似X时有N为有效数字
算法设计的原则
xi1

xi

xi2 25 2xi
取x0=2代入上式，得x1=7.25，继续递推， 依次得5.35、5.0114、5.000001、5.0000000001 …
牛顿法
function x=maqnewt(fun,x0,x1,ep,N) if nargin<5,N=500;end if nargin<4,ep=1e-4;end k=0;
我们可以通过多种方法将方程式
f x 0 转化为 x gx
例如方程 x2 c 0, c 0
可以转化为以下不同形式
(1) x x2 x c
(2)
x c x
(3) x x x2 c 1 x c 2x 2 x
不动点迭代法
从给定的初值x0，按上式可以得到一个数列：
r= 1.4656 -0.2328 + 0.7926i -0.2328 - 0.7926i
>>polyval(c, r(1)) ans = -2.5535e-015
非线性方程求解函数fzero
fzero
对于一般的单个超越方程，可以采用fzero函数求解 fzero函数结合使用二分法、割线法和可逆二次内插
法
函数fzero
[x,fval,exitflag,output] = fzero(fun,x0,options, p1, p2, ...)
此函数的作用求函数fun在x0附件的零值点，x0是 标量
x
数学模型：范德华方程变形可得关于V的非线性 方程
f
(V
)

(p

an2 V2
)(V

nb)

nRT

0
非线性方程
非线性方程包括：高次代数方程、超越方程及其它们 的组合
与线性方程相比，非线性方程求解问题无论从理论上 还是从计算公式上都要复杂得多
对于高次代数方程，当次数>4时，则没有通解公式可 用，对于超越方程既不知有几个根，也没有同样的求 解方式。实际上，对于n≥3代数方程以及超越方程都 采用数值方法求近似根。
非线性方程数值求解原理
逐步扫描法
开始
y
输入初值A，步长h
X1= A，Y1= f(X1)
X2=X 1+h Y2=f (X2)
Y1*Y2 > 0
否 打印A，h，X1
X1= X2 Y1=Y2
是
y = f (x)
a+h
0a
a+2h
y
b
x
y = f (x)
结束
0
a
b
x
逐步扫描法效率较低，常用于求根的初始近似值
不动点迭代的图形解释
y
y
y=x
y=x
(p1,p1)
(p0,g(p0))
P
y = g(x)
P
y = g(x)
(p1,p1)
(p0,g(p0))
O
Pp2 p1
p0
x
0 g 'P 1
单调收敛
O p1 Pp2
p0
x
1 g 'P 0 振荡收敛
不动点迭代的图形解释
y
y = g(x)
y= x
牛顿法注意事项
在有根区间[a,b]上，f '(x) 0, f "(x) 连续且不变号，则只 要选取的初始近似根x0满足 f (x0 ) f "(x0) 0 ，切线法 必定收敛。
在单根附近，牛顿公式恒收敛，而且收敛速度很快。 但是需要注意如果初始值不在根的附近，牛顿公式 不一定收敛
在实际使用中，牛顿法最好与逐步扫描法结合起来， 先通过逐步扫描法求出根的近似值，然后用牛顿公 式求其精确值，以发挥牛顿法收敛速度快的优点
a1=b1; b1=a1+h; k=k+1; end for i=1:k
if x1(i)-x2(i)~=0 [x1(i),x2(i)] end end
求方程根的精确解
非线性方程(组)的求解一般采用迭代法进行。 迭代法是一种重要的逐次逼近方法。这种方法
用某个固定公式反复校正根的近似值，使之逐 步精确化，最后得到满足精度要求的结果。 常见的算法有
while abs(feval(fun,x))>ep|(b-
a>ep)
if feval(fun,x)*feval(fun,a)<0
b=x;
else
a=x;
end
x=(a+b)/2.0;k=k+1;
end
二分法是一种可靠的算法，但计算速度较慢
二分法计算示例
计算方程 x2 2 0 的正数解
不动点迭代法
函数roots
r = roots(c)，用于求解多项式的根 其中，行向量c的元素是多项式的系数，按多
项式次数降序排列 如果c中含有n+1个元素，则多项式为n次 roots可以获得多项式的所有根
例题：
求方程 x3 x2 1 的根
>>c = [1 -1 0 -1]; >>r = roots(c)
例题：f x x2 2x 3 0
构造迭代格式如下：
1）x 2x 3,即x 2x 3，取x0 4,
得迭代格式如下：
xk1 2xk 3，k 0,1,2,
求得x1 2x0 3 11 3.317 x2 2x1 3 9.634 3.104
则 f x对x x0的泰勒级数展开式为
f
x
f
x0 x x0
f
'

x0



x
x0 2!
2
f
" x0
f x f x0 x x0 f ' x0 0
x
x0

f x0 f ' x0
xk 1

xk
寻找根区间 function masearch(fun,a,b,h) n=(b-a)/h; x1=zeros(1,n); x2=zeros(1,n); a1=a; b1=a1+h; k=1;
while(b1<b) if feval(fun,a1)*feval(fun,b1)<0 x1(k)=a1; x2(k)=b1; else a1=b1; b1=a1+h; continue; end
(p0,g(p0))
(p1,p1)
P
O
P p0 p1
p2
x
1 g'P
单调发散
y
y = g(x)
y=x
P
(p0,g(p0))
(p1,p1)
O
p1 P p0 p2
x
g 'P 1
振荡发散
迭代法
function x=maiter(phi,x0,ep,N) if nargin<4 N=500;end if nargin<3 ep=1e-4;end k=0; while k<N x=feval(phi,x0); if abs(x-x0)<ep break; end x0=x;k=k+1; end if k==N, warning(’已达迭代次数上限’); end