MATLAB用二分法、不动点迭代法及Newton迭代（切线）法求非线性方程的根
作者：凯鲁嘎吉 - 博客园
http://www.cnblogs.com/kailugaji/
一、实验原理
二、实验步骤
三、实验过程
1.(程序)
（1）二分法：求 在区间（1，2）之间的根，取
（a）bipart.m:
function [x,m]=bipart(fun,a0,b0,tol)
a=a0;b=b0;
m=1+round(round(log((b-a)/tol))/log(2));
for k=1:m
p=(a+b)/2;
if fun(p)*fun(b)<0
a=p;
else
b=p;
end
x=p;
end
(b)fun1.m:
function f=fun1(x) f=x^3+10*x-20;
（2）不动点迭代法：求方程在附近的根，取
2020/10/17 注：这里求解用的牛顿迭代法，不是不动点迭代法，不动点迭代法的MATLAB程序请参看：MATLAB实例：不动点迭代法求一元函数方程的根 - 凯鲁嘎吉 - 博客园
（a）budong.m:
function [x,k]=budong(fun,x0,tol,m)
for k=1:m
x=fun(x0);
if abs(x-x0)<tol
break;
end
x0=x;
end
x=vpa(x,8);
(b)fun.m
function t=fun(x1) syms x; f=x^3-2*x-5; s=subs(diff(f,x),x,x1); x=x1; f=x^3-2*x-5; t=x-f/s;
（3）牛顿迭代法：求方程在附近的根，取
newton.m:
function x1=newton(t1,esp,m)
syms x;
fun=x^3+2*x-5;
for k=1:m
if abs(subs(diff(fun,\'x\'),x,t1))<esp
x1=t1;
break;
else
if subs(diff(fun,\'x\',2),x,t1)==0
break;
disp(\'解题失败！\')
else
t0=t1;
t1=t0-subs(fun,x,t0)/subs(diff(fun,\'x\'),x,t0);
if abs(t1-t0)<esp
x1=t1;
break;
end
end
end
end
x1=vpa(x1,8);
2.(运算结果)
（1）二分法：
>> [x,m]=bipart(@fun1,1,2,0.0001)
x =
1.5945
m =
14
(2)不动点迭代法：
>> [x,k]=budong(@fun,2,1e-5,100)
x =
2.0945515
k =
4
(3)牛顿迭代法：
>> x1=newton(2,1e-4,20)
x1 =
1.3282689
3.(拓展（方法改进、体会等）)
对于方程的根为重根的情形，newton法求重根只是线性收敛，迭代缓慢，如果对于求重根的情形，对newton法进行改进，取
，
则 。用迭代法
求m重根，则具有二阶收敛性，但要知道的重数m。
计算方程的根是二重根，用newton法与改进方法求根。
源程序：
newton_biroot.m:
function t=newton_biroot(x1) syms x; f=x^4-4*(x^2)+4; s=subs(diff(f,x),x,x1); x=x1; f=x^4-4*(x^2)+4; t=x-f/s;
biroot1.m:
function t=biroot1(x1) syms x; f=x^4-4*(x^2)+4; s=subs(diff(f,x),x,x1); x=x1; f=x^4-4*(x^2)+4; t=x-2*f/s;
budong.m:
function [x,k]=budong(fun,x0,tol,m)
for k=1:m
x=fun(x0);
if abs(x-x0)<tol
break;
end
x0=x;
x=vpa(x,8)
end
x=vpa(x,8);
运行结果：取初值为2
|
k |
xk |
newton法 |
改进方法 |
|
1 |
x1 |
1.75 |
1.5 |
|
2 |
x2 |
1.5982143
|
1.4166667
|
|
3 |
x3 |
1.5115099
|
1.4142157
|
|
4 |
x4 |
1.4644275
|
1.4142157
|
计算4步，改进方法就已经收敛，而newton法只是线性收敛，要达到同样精度需迭代17次。
附结果：
>> [x,k]=budong(@biroot1,2,1e-5,3)
x =
1.5
x =
1.4166667
x =
1.4142157
x =
1.4142157
k =
3
>> [x,k]=budong(@biroot1,2,1e-5,10)
x =
1.5
x =
1.4166667
x =
1.4142157
x =
1.4142136
k =
4
>> [x,k]=budong(@newton_biroot,2,1e-5,50)
x =
1.75
x =
1.5982143
x =
1.5115099
x =
1.4644275
x =
1.439751
x =
1.4270955
x =
1.4206836
x =
1.4174559
x =
1.4158366
x =
1.4150256
x =
1.4146197
x =
1.4144166
x =
1.4143151
x =
1.4142643
x =
1.414239
x =
1.4142263
x =
1.4142199
k =
17
