1. 拉格朗日插值

%Lagrange.m
function y=Lagrange(x,f,x0) 
% 给定的一系列数据点x,f
% x0是我们要预测的值，由于可以有多个，因此用向量表示
% y返回我们的估计值，由于可以有多个，因此用向量表示
        n=length(x); %得出数据点的个数
        y=zeros(n); %初始化,并赋初值0
        for k=1:n %外循环k
            l=1.0 %置l(x)的初值
            for j=1:n %内循环j
                if j~=k %除去j=k的点
                    l=l*(x0-x(j))/(x(k)-x(j)); %累乘得l(x)
                end
            end
            y=y+l*f(k);%累加得y
        end
end

2. 牛顿插值

%Newton.m
function y=Newton(x,xi,yi)
%xi: 离散样点的横坐标值
%yi: 离散样点的纵坐标值
%x: 插值多项式中自变量符号
%y: Newton插值多项式
        n=length(x); %得出数据点的个数
        f=zeros(n,n); %初始化,并赋初值0
        %对差商表的第一列赋值
        for k=1:n
            f(k)=yi(k);
        end
        %求差商表
        for i=2:n %差商表从0阶开始；但是矩阵是从1维开始存储
            for k=1:n
                f(k,i)=(f(k,i-1)-f(k-1.i-1))/(xi(k)-xi(k+1-i));
            end
        end
        %求插值多项式
        y=0;
        for k=2:n
            t=1;
            for j=1:k-1
                t=t*(x-xi(j));
            end
            y=y+f(k,k)*t;
        end
        y=f(1,1)+y;
end

3. 改进欧拉法

%Euler.m
function []=Euler(h,x0,y0,N)
    for i=1:1/h
        x=x0+(i-1)*h;
        K1=f(x0,y0);
        K2=f(x+h,y+h*K1);
        y=y+(h/2)*(K1+K2);
    end
end

4. 经典四阶龙格-库塔方法

%R_K.m
function []=R_K(h,x0,y0,N)
    for i=1:1/h
        x=x0+(i-1)*h;
        K1=f(x0,y0);
        K2=f(x+h/2,y+(h/2)*K1);
        K3=f(x+h/2,y+(h/2)*K2);
        K4=f(x+h,y+h*K3);
        y=y+(h/6)*(K1+2*K2+2*K3+K4);
    end
end

5. 牛顿迭代法

%N_Iteration
function x=N_Iteration(f,df,x0,eps,N)
    k=1;
    while k<N
        if feval(df,x0)==0
            I=-1;
            break;
        else
            x=x0-feval(f,x0)/feval(df,x0);
        end
        if abs(x-x0)<eps 
            I=0;
            break;
        end
        x0=x;
        k=k+1;
    end
end

6. 高斯-塞德尔迭代

%Gauss_Seidel.m
function Gauss_Seidel(A,b,n,x0,tol,N)
    x=zeros(n,1);  % 给x赋值
    k=1;
    while k<N
        for i=1:n
            if i==1
                x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
            elseif i==n
                x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);
            else
                x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i)       
            end  
        end
        if norm(x-x0)<tol
            break;
        end
        x0=x;
        k=k+1;
    end
end

7. 列主元高斯消去法

function x=gauss_column(A,b)      %输入矩阵A和列向量b，返回解向量x
    [ni,nj]=size(b);
    if rank(A)~=rank([A,b])       %若系数矩阵秩和增广矩阵秩不相等，则无解
        fprintf('无解\n')
        return
    else if rank(A)<ni            %若系数矩阵秩和增广矩阵秩相等，但是其秩小于未知量个数，则无穷解
        fprintf('无穷解\n')
        return
    else        
        for j=1:ni
            [tv,ti]=max(A(j:ni,j));             %找出该列中按模最大的元素
            A([ti+j-1,j],:)=A([j,ti+j-1],:);    %交换行
            for i=j+1:ni                        %消去过程
                d=-A(i,j)/A(j,j);
                A(i,:)=A(i,:)+d*A(j,:);
                b(i)=b(i)+d*b(j);
            end
        end
        x=zeros(size(b));                       %初始化解向量
        x(ni)=b(ni)/A(ni,ni);
        for i=ni-1:-1:1                        %回代过程
            x(i)=(b(i)-sum(x.*A(i,:)'))/A(i,i);
        end
    end
end