Matlab(二分法、简单迭代法、牛顿迭代法、弦截法)&&(高斯消去法、矩阵的三角分解法、雅可比迭代法、高斯-赛德尔迭代法)

迈不过友情╰ 2022-04-24 18:48 474阅读 0赞

function root= Secant( f,a,b,eps )
    if(nargin==3)
        eps=1.0e-4;
    end
    f1=subs(sym(f),findsym(sym(f)),a);
    f2=subs(sym(f),findsym(sym(f)),b);
    if(f1==0)
        root=a;
    end
    if(f2==0)
        root=b;
    end
    if(f1*f2>0)
        disp('两端点乘积大于零');
        return;
    else
        tol=1;
        fa=subs(sym(f),findsym(sym(f)),a);
        fb=subs(sym(f),findsym(sym(f)),b);
        root=a-(b-a)*fa/(fb-fa);
        while(tol>eps)
            r1=root;
            fx=subs(sym(f),findsym(sym(f)),r1);
            s=fx*fa;
            if(s==0)
                root=r1;
            else
                if(s>0)
                    root=b-(r1-b)*fb/(fx-fb);
                else
                    root=a-(r1-a)*fa/(fx-fa);
                end
            end
            tol=abs(root-r1);
        end
    end
            
    
    end
    
    ------------------------------------------------------
    function [ p1,err,k,y ] = Newtonroot1( f,df,p0,eps,maxl )
    p0,feval('f',p0)
    for k=1:maxl
        p1=p0-feval('f',p0)/feval('df',p0);
        err=abs(p1-p0);
        p0=p1;
        p1,err,k,y=feval('f',p1)
        if(err<eps)|(y==0)
            break;
        end
    
    end
    
    ------------------------------------------------------------
    function [ p0,k,err,p ] = Fixpt( g,p0,tol,maxl )
    P(1)=p0;
    for k=2:maxl    
    P(k)=feval('g',P(k-1));
    k,err=abs(P(k)-P(k-1))
    p=P(k);
    if(err<tol)
        break;
    end;
    if(k==maxl)
        disp('over');
    end
    
    
    end
    P
    
    function  root  = HalfInterval( f,a,b,eps )
    if(nargin==3)
        eps=1.0e-4;
    end
    f1=subs(sym(f),findsym(sym(f)),a);
    f2=subs(sym(f),findsym(sym(f)),b);
    if(f1==0)
        root=a;
    end
    if(f2==0)
        root=b;
    end
    if(f1*f2>0)
        disp('两端点函数值乘积大于0！');
        return;
    else
        root=FindRoots(f,a,b,eps);
    end
      
                        
                   
    
    end
     
    function r=FindRoots(f,a,b,eps)
            f3=subs(sym(f),findsym(sym(f)),a);
             f4=subs(sym(f),findsym(sym(f)),b);
            mf=subs(sym(f),findsym(sym(f)),(a+b)/2);
            if(f3*mf>0)
                t=(a+b)/2;
                r=FindRoots(f,t,b,eps);
            else
                if(f3*mf==0)
                    r=(a+b)/2;
                else
                    if(abs(b-a)<=eps)
                        r=(b+3*a)/4;
                    else
                        s=(a+b)/2;
                        r=FindRoots(f,a,s,eps);
                    end
                end
            end
    -----------------------------------------------------------------
           function [ x,det,flag ] = Guass( A,b )
    %UNTITLED3 Summary of this function goes here
    %   Detailed explanation goes here
    [n,m]=size(A);
    nb=length(b);
    if(n~=m)
        error('the rows and columns of matrix A must equal');
        return;
    end
    if(m~=nb)
        error('the rows and columns of  A must equal b');
        return;
    end
    flag='OK';
    det=1;
    x=zeros(n,1);
    for(k=1:n-1)
        max1=0;
        for(i=k:n)
        if(abs(A(i,k))>max1)
            max1=abs(A(i,k));
            r=i;
        end
       end
    if(max1<1e-10)
        flag='failure';
        return;
    end
    if(r>k)
        for(j=k:n)
            z=A(k,j);
            A(k,j)=A(r,j);
            A(r,j)=z;
        end
        z=b(k);
        b(k)=b(r);
        b(r)=z;
        det=-det;
    end
    for(i=k+1:n)
        m=A(i,k)/A(k,k);
        for(j=k+1:n)
            A(i,j)=A(i,j)-m*A(k,j);
        end
        b(i)=b(i)-m*b(k);
    end
    det=det*A(k,k);
    end
    det=det*A(n,n);
    if(abs(A(n,n))<1e-10)
        flag='failure';
        return;
    end
    for(k=n:-1:1)
        for(j=k+1:n)
            b(k)=b(k)-A(k,j)*x(j);
        end
        x(k)=b(k)/A(k,k);
    end
    x(k)=b(k)/A(k,k);
    end
    
    --------------------------------------
    function [ L,U,flag ] = LU_Demo( A )
    %UNTITLED Summary of this function goes here
    %   Detailed explanation goes here
    [n,m]=size(A);
    if(n~=m)
    return;
    end
    L=eye(n);
    U=zeros(n);
    flag='ok';
    for(k=1:n)
        for(j=k:n)
            z=0;
            for(q=1:k-1)
                z=z+L(k,q)*U(q,j);
            end
            U(k,j)=A(k,j)-z;
        end
        if(abs(U(k,k))<eps)
            flag='f';
            return;
        end
        for(i=k+1:k-1)
            z=0;
            for(q=1:k-1)
                z=z+L(i,q)*U(q,k);
            end
            L(i,k)=(A(i,k)-z)/U(k,k);
        end
        
    end
    -----------------------------------------------------------
    function [ x,n ] = Jacobi( A,b,x0,eps,varargin )
    %UNTITLED3 Summary of this function goes here
    %   Detailed explanation goes here
    if(nargin==3)
        eps=1.0e-6;
        M=200;
    elseif(nargin<3)
        error
        return
    elseif(nargin==5)
        M=varargin{1};
    end
    D=diag(diag(A));
    L=-tril(A,-1);
    U=-triu(A,1);
    B=D\(L+U);
    f=D\b;
    x=B*x0+f;
    n=1;
    while(norm(x-x0)>=eps)
        x0=x;
        x=B*x0+f;
        n=n+1;
    if(n>=M)
        disp('duo');
        return;
    end
    end
    
    -----------------------------------------------
    function [ x,n ] = GaussSeidel( A,b,x0,eps,M )
    %UNTITLED4 Summary of this function goes here
    %   Detailed explanation goes here
    if(nargin==3)
        eps=1.0e-6;
        M=200;
    elseif(nargin==4)
        M=200;
    elseif(nargin<3)
        error;
        return;
    end
    D=diag(diag(A));
    L=-tril(A,-1);
    U=-triu(A,1);
    G=(D-L)\U;
    f=(D-L)\b;
    x=G*x0+f;
    n=1;
    while(norm(x-x0)>=eps)
        x0=x;
        x=G*x0+f;
        n=n+1;
    if(n>=M)
        disp('duo');
        return;
    end
    end