数值分析算法 MATLAB 实践 非线性方程(组)求解
4.1 二分法
function [c,err,yc]=bisectmethod(fun,a,b,delta)
% 二分法function [c,err,yc]= bisectmethod (fun,a,b,delta)
%a为左区间，b为右区间，delta为区间误差限
%c求解的值，err误差，yc:c在f(x)处的值
ya=feval('fun',a);
yb=feval('fun',b);
if yb==0
c=b;
return;
end
if(ya*yb)>0
disp('(a,b)不是有根区间');
return;
end
max1=1+round((log(b-a)-log(delta))/log(2));
for k=1:max1
c=(a+b)/2;
yc=feval('fun',c);
if yc==0
a=c;
b=c;
return;
elseif(yb*yc)>0
b=c;
yb=yc;
else
a=c;
ya=yc;
end
if(b-a)<delta
break;
end
end
c=(a+b)/2;
err=abs(b-a);
yc=feval('fun',c);
end
% 二分法function [c,err,yc]=bisectmethod(fun,a,b,delta)
%a为左区间，b为右区间，delta为区间误差限
%c求解的值，err误差，yc:c在f(x)处的值
eps = 1e-4;
[x0,err,y_x0]=bisectmethod(@fun,-2,-1,eps);
disp(['方程组的解：x0= ']);
disp(x0)
function [xc,k_cnt] = bisect(f,a,b,tol)
% 二分法function [xc,k_cnt] = bisect(f,a,b,tol)
%a为左区间，b为右区间，tol为区间误差限
%xc求解的值
if sign(f(a))*sign(f(b))>=0
error('f(a)f(b)<0 not satisfied!') % 停止运行
end
fa = f(a);
fb = f(b);
k_cnt=0;
while (b-a)/2 > tol
c = (a+b)/2;
fc = f(c);
if fc == 0
break
end
if sign(fc)*sign(fa) < 0
% [a,c]是新的区间
b = c;
fb = fc;
else % [b,c]是新的区间
a = c;
fa = fc;
end
k_cnt=k_cnt+1;
end
xc = (a+b)/2;
% 新的中点就是最优估计
end
% 二分法function [xc,k_cnt] = bisect(f,a,b,tol)
%a为左区间，b为右区间，tol为区间误差限
%xc求解的值
[x1,k1_cnt] = bisect(@fun,-2,-1,eps);
disp(['迭代次数：k1_cnt= ']);
disp(k1_cnt)
disp(['方程组的解：x1= ']);
disp(x1)
function y=fun(x)
y=x.^3-3*x-1;
end
4.2 牛顿迭代法
function[p1,err,k,y]=NewtonFunmethod(fun,dfun,p0,delta,max1)
%fun是给定的非线性函数 %p0为初始值 %delta为给定误差界
%max1迭代次数的上限 %p1为所求得的方程的近似解 %err为p1-p0的绝对值
%k为所需要的迭代次数 %y=fun(p1)
for k=1:max1
p1 = p0-fun(p0)/dfun(p0);
%一定要有fun,dfun .m文件
%p1=p0feval('fun',p0)/feval('dfun',p0);
err=abs(p1-p0);
p0=p1;
%y=feval('fun',p1);
y= fun(p1);
if(err<delta)|(y==0)
break;
end
end
end
%function [p1,err,k,y]=NewtonFunmethod(fun,dfun,p0,delta,max1)
%NewtonFunmethod对于初值比较敏感，选择不好，出现不收敛情况
delta = 1e-7;
p0 = 1.2;
it_max = 100;
[x2,err,k2_cnt,y]=NewtonFunmethod(@fun1,@dfun1,p0,delta,it_max);
disp(['迭代次数：k2_cnt= ']);
disp(k2_cnt)
disp(['方程组的解：x2= ']);
disp(x2)
function y=fun1(x)
y=x.^3-3*x+2;
end
function y=dfun1(x)
y=3*x.^2-3;
end
4.3 割线法
function[p1,err,k,y]=secantFunmethod(fun,p0,p1,delta,max1)
%割线法fun是给定的非线性函数 %p0，p1为初始值 %delta为给定误差界
%max1迭代次数的上限 %p1为所求得的方程的近似解 %err为p1-p0的绝对值
%k为所需要的迭代次数 %y=fun(p1)
for k=1:max1
p2 = p1-fun(p1)*(p1-p0)/(fun(p1)-fun(p0));
err=abs(p2-p1);
p0=p1;
p1=p2;
y = fun(p1);
if(err<delta)||(abs(y)<0.00001)
break;
end
end
delta = 1e-6;
p0 = 1.0;
p1 = 1.5;
it_max = 100;
[x3,err,k3_cnt,y]=secantFunmethod(@fun2,p0,p1,delta,it_max);
disp(['迭代次数：k3_cnt= ']);
disp(k3_cnt)
disp(['方程组的解：x3= ']);
disp(x3)
function y=fun2(x)
y=x.^3+4*x.^2-10;
end
4.4 非线性方程组-牛顿法
%牛顿迭代法 计算非线性方程组
%输入 x0 为迭代初值
%tol 为误差容限 如果缺省 默认为 10 的-10 次方
%data 用来存放计算的中间数据便于计算收敛情况分析
function [x,n,data]=new_ton(x0,tol)
if nargin==1
tol=1e-10;
end
x1=x0-f1(x0)/df1(x0);
n=1;
%迭代过程
while (norm(x1-x0)>tol)
x0=x1;
x1=x0-f1(x0)/df1(x0);
n=n+1;
%data 用来存放中间数据
data(:,n)=x1;
end
x=x1;
end
%牛顿迭代法的 方程函数
function f=f1(x0)
x=x0(1);
y=x0(2);
f1=x^2-2*x-y+0.5;
f2=x^2+4*y^2-4;
%最后方程函数 以行向量输出
f=[f1 f2];
end
% 求解非线性方程表达式的雅可比矩阵
function f=df1(x0)
x=x0(1);
y=x0(2);
f=[2*x-2 -1
2*x 8*y];
end
x1=[1 1];
[x_1,k1_cnt,data]=new_ton(x1);
disp(['迭代次数：k1_cnt= ']);
disp(k1_cnt)
disp(['方程组的解：x_1= ']);
disp(x_1)
%抽取 data1 中第一个变量数据 画出曲线
subplot(2,1,1)
plot(data(1,:)),title('x 在迭代中的变化')
%抽取 data 中的第二个变量数据 画出其变化曲线
subplot(2,1,2)
plot(data(2,:)),title('y 在迭代中的变化')
function [allx,ally,r,n]=MultiNewton(F,x0,eps)
% allx是用来记录每一步迭代的点的矩阵，ally是用来记录迭代点对应计算出来的函数值的，
% r是满足精度要求的最终迭代点 n迭代次数
if nargin==2
eps=1.0e-4;
end
x0 = transpose(x0);
Fx = subs(F,transpose(symvar(F)),x0);
var = transpose(symvar(F));
dF = jacobian(F,var);
dFx = subs(dF,transpose(symvar(F)),x0);
n=dFx;
r=x0-inv(dFx)*Fx';
n=1;
tol=1;
N=100;
symx=length(x0);
ally=zeros(symx,N);
allx=zeros(symx,N);
while tol>eps
x0=r;
Fx = subs(F,transpose(symvar(F)),x0);
dFx = subs(dF,transpose(symvar(F)),x0);
r=vpa(x0-inv(dFx)*Fx');
tol=norm(r-x0)
if(n>N)
disp('迭代步数太多，可能不收敛！');break;
end
allx(:,n)=x0;
ally(:,n)=Fx;
n=n+1;
end
function f = Multifun(x)
% 只需要将fun.m文件中的 k更改为所需维数，文件中的方程替换为你自己所需的方程
% x(1)->x x(2)->y x(3)->z
%{
f(1)=3*x(1)-cos(x(1)*x(2))-1/2;
f(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
f(3)= exp(-1*x(1)*x(2))+20*x(3)+(10*pi-3)/3;
%}
k=3;
for i=1:k
x(i)=sym (['x',num2str(i)]);
end
f(1)=3*x(1)-cos(x(1)*x(2))-1/2;
f(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
f(3)= exp(-1*x(1)*x(2))+20*x(3)+(10*pi-3)/3;
end
allx0 = [1,1,1];
[allx,ally,r,k_cnt]=MultiNewton(Multifun,allx0);
disp(['迭代次数：k_cnt= ']);
disp(k_cnt)
disp(['方程组的解：allx= ']);
disp(allx)
