math_exp41.m

matlab 牛顿法迭代,为数值分析课上的作业,用借鉴于其他的迭代计算中.

clear
clc
%%%   matlab 提供的文件 %%%%%%%%%%%%%%%%%%%%%%%%% 
x1=[1 1 1]';
x2=[0 0 0]';
%x1=[100 100 100]';
%x1=[20 20 10]';
%x2=x1;
maxcycle=1000;
tic
[Xf,Fval,ExitFlag,Iters]=mmfsolve('fun1',x1); % Broyden
time11=toc;
tic
[Xf1,Fval1,ExitFlag1,Iters1]=mmfsolve('fun2',x2); % Broyden
time12=toc;

tol=[1e-7, 1e-7]; params=[1000, 1, 0]; % Newton
tic
[sol, it_hist, ierr] = nsol(x1, 'fun1', tol, params); % Newton
time01=toc;
tic
[sol1, it_hist1, ierr1] = nsol(x2, 'fun2', tol, params); % Newton
time02=toc;

%0000000000000000 自己编程写特定方程的牛顿法和Broyden法
e=10^-7;
tic
for cycle=1:maxcycle
xx=x1;
% -----------------  Newton     ---------------------
%------------------  function 1 ---------------------
k=1;
while (k<100)
    ff1=12*xx(1,k)-xx(2,k)^2-4*xx(3,k)-7;
    ff2=xx(1,k)^2+10*xx(2,k)-xx(3,k)-11;
    ff3=xx(2,k)^3+10*xx(3,k)-8;  
    FF=[ff1; ff2; ff3];
    dFF=[12 -2*xx(2,k) -4; 2*xx(1,k) 10 -1; 0 3*xx(2,k)^2 10];
    dxx=inv(dFF)*FF;
    xx(:,k+1)=xx(:,k)-dxx;
    if norm(dxx,1) < e
        break;
    end
     k=k+1;
end
iterative1=k;
end
time1=toc;

tic
for cycle=1:maxcycle
x=x2;
%------------------- function 2-----------------------
k=1;
while (k<100)    
    f1=3*x(1,k)-cos(x(2,k)*x(3,k))-0.5;    
    f2=x(1,k)^2-81*(x(2,k)+0.1)^2+sin(x(3,k))+1.06;    
    f3=exp(-x(1,k)*x(2,k))+20*x(3,k)+(10*pi-3)/3;
    F=[f1; f2; f3];
    dF=[3 sin(x(3,k)) sin(x(2,k)); 2*x(1,k) -162*(x(2,k)+0.1) cos(x(3,k)); 
        -x(2,k)*exp(-x(1,k)*x(2,k)) -x(1,k)*exp(-x(1,k)*x(2,k)) 20];
    dx=inv(dF)*F;
    x(:,k+1)=x(:,k)-dx;
    if norm(dx,1)<e
        break;
    end;
    k=k+1;
end;
iterative2=k;
end
time2=toc;
% -----------------  Broyden秩1法和   逆Broyden秩1方法  ---------------------
%------------------  function 1 ---------------------
tic
for cycle=1:maxcycle
bxx=x1;
%H=inv([12 -2*bxx(2,1) -4; 2*bxx(1,1) 10 -1; 0 3*bxx(2,1)^2 10]);%初始化H0
A=eye(3);
for k=1:100
bf1=12*bxx(1,k)-bxx(2,k)^2-4*bxx(3,k)-7;
bf2=bxx(1,k)^2+10*bxx(2,k)-bxx(3,k)-11;
bf3=bxx(2,k)^3+10*bxx(3,k)-8;
F(:,k)=[bf1; bf2; bf3];
bxx(:,k+1)=bxx(:,k)-inv(A)*F(:,k);
%bxx(:,k+1)=bxx(:,k)-H*F(:,k);
bf11=12*bxx(1,k+1)-bxx(2,k+1)^2-4*bxx(3,k+1)-7; 
bf21=bxx(1,k+1)^2+10*bxx(2,k+1)-bxx(3,k+1)-11;
bf31=bxx(2,k+1)^3+10*bxx(3,k+1)-8;
F(:,k+1)=[bf11; bf21; bf31];
    y=F(:,k+1)-F(:,k);
    s1=bxx(:,k+1)-bxx(:,k);
    dbxx=s1;
     A=A+(y-A*s1)*s1'/(s1'*s1); 
%   H=H+(s1-H*y)*s1'*H/(s1'*H*y); 
    if norm(dbxx,1)<e
        break;
    end;
end;
iterative3=k;
end
time3=toc;

tic
for cycle=1:maxcycle
bxx2=x1;
H=inv([12 -2*bxx2(2,1) -4; 2*bxx2(1,1) 10 -1; 0 3*bxx2(2,1)^2 10]);%初始化H0
%A=eye(3);
for k=1:100
bf1=12*bxx2(1,k)-bxx2(2,k)^2-4*bxx2(3,k)-7;
bf2=bxx2(1,k)^2+10*bxx2(2,k)-bxx2(3,k)-11;
bf3=bxx2(2,k)^3+10*bxx2(3,k)-8;
F(:,k)=[bf1; bf2; bf3];
%bxx(:,k+1)=bxx(:,k)-inv(A)*F(:,k);
bxx2(:,k+1)=bxx2(:,k)-H*F(:,k);
bf11=12*bxx2(1,k+1)-bxx2(2,k+1)^2-4*bxx2(3,k+1)-7; 
bf21=bxx2(1,k+1)^2+10*bxx2(2,k+1)-bxx2(3,k+1)-11;
bf31=bxx2(2,k+1)^3+10*bxx2(3,k+1)-8;
F(:,k+1)=[bf11; bf21; bf31];
    y=F(:,k+1)-F(:,k);
    s1=bxx2(:,k+1)-bxx2(:,k);
    dbxx=s1;
%      A=A+(y-A*s1)*s1'/(s1'*s1); 
  H=H+(s1-H*y)*s1'*H/(s1'*H*y); 
    if norm(dbxx,1)<e
        break;
    end;
end;
iterative5=k;
end
time32=toc;

%--------------------- function 2 ---------------
tic
for cycle=1:maxcycle
bx=x2;
%H=inv([3 sin(bx(3,1)) sin(bx(2,1));
%    2*bx(1,1) -162*(bx(2,1)+0.1) cos(bx(3,1));
%    -bx(2,1)*exp(-bx(1,1)*bx(2,1)) -bx(1,1)*exp(-bx(1,1)*bx(2,1)) 20]);

A=eye(3);
for k=1:100
    f1=3*bx(1,k)-cos(bx(2,k)*bx(3,k))-0.5;
    f2=bx(1,k)^2-81*(bx(2,k)+0.1)^2+sin(bx(3,k))+1.06;
    f3=exp(-bx(1,k)*bx(2,k))+20*bx(3,k)+(10*pi-3)/3;
    F(:,k)=[f1; f2; f3];
    bx(:,k+1)=bx(:,k)-inv(A)*F(:,k);
    %bx(:,k+1)=bx(:,k)-H*F(:,k);
    f11=3*bx(1,k+1)-cos(bx(2,k+1)*bx(3,k+1))-0.5;
    f21=bx(1,k+1)^2-81*(bx(2,k+1)+0.1)^2+sin(bx(3,k+1))+1.06;
    f31=exp(-bx(1,k+1)*bx(2,k+1))+20*bx(3,k+1)+(10*pi-3)/3;
    F(:,k+1)=[f11; f21; f31];
    y=F(:,k+1)-F(:,k);
    s=bx(:,k+1)-bx(:,k);
    dbx=s;
   A=A+(y-A*s)*s'/(s'*s); 
    %H=H+(s-H*y)*s'*H/(s'*H*y);
    if norm(dbx,1)<e
        break;
    end;
end;
iterative4=k;
end
time4=toc;

tic
for cycle=1:maxcycle
bx2=x2;
H=inv([3 sin(bx2(3,1)) sin(bx2(2,1));
   2*bx2(1,1) -162*(bx2(2,1)+0.1) cos(bx2(3,1));
   -bx2(2,1)*exp(-bx2(1,1)*bx2(2,1)) -bx2(1,1)*exp(-bx2(1,1)*bx2(2,1)) 20]);

%A=eye(3);
for k=1:100
    f1=3*bx2(1,k)-cos(bx2(2,k)*bx2(3,k))-0.5;
    f2=bx2(1,k)^2-81*(bx2(2,k)+0.1)^2+sin(bx2(3,k))+1.06;
    f3=exp(-bx2(1,k)*bx2(2,k))+20*bx2(3,k)+(10*pi-3)/3;
    F(:,k)=[f1; f2; f3];
    %bx(:,k+1)=bx(:,k)-inv(A)*F(:,k);
    bx2(:,k+1)=bx2(:,k)-H*F(:,k);
    f11=3*bx2(1,k+1)-cos(bx2(2,k+1)*bx2(3,k+1))-0.5;
    f21=bx2(1,k+1)^2-81*(bx2(2,k+1)+0.1)^2+sin(bx2(3,k+1))+1.06;
    f31=exp(-bx2(1,k+1)*bx2(2,k+1))+20*bx2(3,k+1)+(10*pi-3)/3;
    F(:,k+1)=[f11; f21; f31];
    y=F(:,k+1)-F(:,k);
    s=bx2(:,k+1)-bx2(:,k);
    dbx=s;
   %A=A+(y-A*s)*s'/(s'*s); 
    H=H+(s-H*y)*s'*H/(s'*H*y);
    if norm(dbx,1)<e
        break;
    end;
end;
iterative6=k;
end
time42=toc;