MATLAB求解⾮线性⽅程组的五种⽅法
MATLAB求解⾮线性⽅程组的五种⽅法
求解线性⽅程分为两种⽅法–⼆分法和迭代法
常见的⽅法⼀共有5种
⼆分法
迭代法
⽜顿法
割线法
拟⽜顿法
Halley法
使⽤条件
⼆分法需要知道两个⾃变量，分别是⼀个根的两侧
⽜顿法迭代法是最常⽤的⽅法，收敛性信赖于初值，取不同的初值可以的⽅程不同的根，函数⽤的是⼀阶导数，输⼊的是⼀个猜想的可能的值
割线法给定两个初值再带⼊计算,⽐如要在2附近求⼀个根，那就可以假设这个范围是（1.9，2）
拟⽜顿法这个⽐较⽅便，⽤时最好可以到⼀个好的初始值
Halley法需要知道函数值以及它的⼀阶求导、⼆阶求导
这⾥我就从计算代码的⾓度来讲解，在下⾯也会按照上⾯这个顺序给出代码，遇到⽅程组直接带⼊已知条件就可以得到答案。
⼆分法
基本函数是这样⼦的：y = dichotomy(fun,a,b,tol)；⼆分法的算法要输⼊四个变量，fun,a,b,tol：函数，⼀个根的左右点，tol=1.0e-6
function y =fun(x)
y = x^3-5* x +4.272;
上⾯这个就是定义的fun，每次的输⼊的⽅程不同，第⼀条不动，直接改第⼆⾏就可以的。⽐如这⾥我们要计算的⽅程y = x^3 - 5 * x + 4.272;
我们是可以通过简单计算得到⼀个根的两侧分别是1和1.3
那在窗⼝指令指令中输⼊x=dichotomy(’fun‘,1,1.3,1.0e-6)就可以得到结果
function y =dichotomy(fun,a,b,tol)
if nargin <4
tol =1.0e-5;
end
n =1;
if feval(fun,a)*feval(fun,b)<0
c =(a+b)/2;
while(abs(b-c)>tol)&&(abs(feval(fun,c))>tol)
if(feval(fun,c)*feval(fun,a)>0)
a = c;
c =(a+b)/2;
elseif(feval(fun,c)*feval(fun,a)<0)
b = c;
c =(a+b)/2;
else
y = c;
tol =100;
end
n = n +1;
end
y = c;
elseif feval(fun,a)==0
y = a;
elseif feval(fun,b)==0
y = b;
else
disp('there may not be a root in the interval');
end
n
function y =fun(x)
y = x^3-5* x +4.272;
⽜顿法
还是⽤刚才那道题，y = x^3 - 5 * x + 4.272，⼀阶导是y = 3 * x^2 - 5;
function y =dfun(x)
y =3* x^2-5;
下⾯的是具体的算法，根据x = newton(x0,tol)，我们只需要输⼊⼀个我们猜想的值就可以。但是有⼀定的误差
function x =newton(x0,tol)
if nargin <2
tol =1.0e-5;
end
x = x0 -fun(x0)/dfun(x0);
n =1;
while(norm(x-x0)>tol)&&(n<1000)
x0 = x;
x = x0 -fun(x0)/dfun(x0);
n = n +1;
end
n
割线法
这⾥我们⽤割线法求y = x^3 - 5 * x + 4.272在⽅程x=2的根，输⼊上要⽤两个初始值，⽐如说现在来计算就可以输⼊x=secant(2,1.9，10e-6)
x = x0 -fun(x0)*(x0 - x1)/(fun(x0)-fun(x1));
n =1;
while(abs(x0-x1)> tol)&&(n <=1000)
x1 = x0;
x0 = x;
x = x0 -fun(x0)*(x0 - x1)/(fun(x0)-fun(x1));
n = n +1;
end
n
拟⽜顿法
这⾥我们可以直接到⼀个初始值输⼊，⽐如说broyden2,10e-6)，如果不知道不确定也没关系，⾄少要知道⼀个范围。⽐如说给个范围（0.5，0.5）有下⾯这个函数
function y =funm(x)
y(1,1)=x(1,1)-0.7*sin(x(1,1))-0.2*cos(x(2,1));
y(2,1)=x(2,1)-0.7*cos(x(1,1))+0.2*sin(x(2,1));
那就可以输⼊x = broyden(x0,tol)
function x =broyden([0.5，0.5],tol)
if nargin <2
tol =1.0e-5;
endmatlab求导
A=eye(size(x0,1));
x = x0 -A \ funm(x0);
n =1;
while(norm(x - x0)> tol)&&(n <1000)
x0 = x;
x = x0 -A \ funm(x0);
p = x - x0;
q =funm(x)-funm(x0);
A=A+(q -A*p)*p'/norm(p)^2;
n = n +1;
end
n
Halley法
这个要求⼆阶导，⽐如说第⼀个道题，y = x^3 - 5 * x + 4.272;，⼆阶导数是下⾯这个输⼊
function y =d2fun(x)
y =6* x;
最后输⼊x = halley(1，10e-6)就可以计算出⼀个结果
m =size(x0,1);
x = x0 -(eye(m)-1/2*(dfun(x0) \ d2fun(x0))*(dfun(x0) \ fun(x0))) \ ... (dfun(x0) \ fun(x0));
n =1;
while(norm(x - x0)> tol)&&(n <1000)
x0 = x;
x = x0 -(eye(m)-1/2*(dfun(x0) \ d2fun(x0))*(dfun(x0) \ fun(x0))) \ ... (dfun(x0) \ fun(x0));
n = n +1;
end
n