MATLAB当中一些简单的数值分析函数总结

求解函数的近似根：
二分法相关的MATLAB函数

function [x]=dichotomy(f,a,b,TOL) 
% f：原函数，a，b：求根区间，需要自定义给出，TOL：误差范围
N=1+fix(log2((b-a)/TOL));
%fix：取整数，并根据|X_N-X*|<=(b-a)/2^N<TOL,求N
n=1;
if (f(a)*f(b)>0)
    disp('根不在输入的区间里，请重新输入区间');
    return
else
    while n <= N
        x=(a+b)/2;
        if f(a)*f(x)>0
            a=x;
        else
            b=x;
        end
        n=n+1;
    end
    disp(x);
end
end

不动点迭代法相关的MATLAB函数
p.s. 请注意，不动点迭代法的原理当中，首先需要对原函数f(x)进行变化，再代入到下列函数中

function [x] = fixPoint(f,x0,TOL)
% f：原函数，x0：初始的迭代值，TOL：误差范围
x1 = f(x0);
while abs(x1 - x0)>TOL
    x0 = x1;
    x1 = f(x0);
end
x = x1;
disp(x)
end

Newton法相关的MATLAB函数

function [X]=Newton(f,p0,TOL,N)
format long;
n=1;%初始迭代次数
syms x;
while n<=N
    if abs(subs(diff(f(x)),x,p0))<TOL
        X=p0;
        break;
    else
        if subs(diff(f(x),2),x,p0)==0
            disp('Method failed');
            break;
        else
            p=p0-f(p0)/subs(diff(f(x)),x,p0);
            p=eval(p);%将exp的值转为小数值
            if(abs(p-p0)<TOL)
                X=p;
                break;
            else
                p0=p;
            end
        end
    end
    n=n+1;
end
disp(X)
end

上述三种函数（二分法，不动点迭代法，牛顿法）主要用于估算函数的根。

求解方程组的近似根：
高斯消去法相关的MATLAB函数

function [x] = gauss (A,b)
%A：系数矩阵，b：等式右侧常数
n = length (b);
for k = 1: n-1
    index = (k+1:n);
    m = -A(index,k)/A(k,k);
    A(index,index) = A(index,index) + m * A(k,index);
    b(index) = b(index) + m *b(k);
end
x = zeros(n,1);
x(n) = b(n)/A(n,n);
for i = n-1:-1:1
    x(i) = ( b(i) - A(i,(i+1:n))*x((i+1:n)))/A(i,i);
end
end

Jacobi迭代法相关的MATLAB函数

function [x] = jacobi(A,b,TOL)
%A：系数矩阵，b：等式右侧常数
D=diag(diag(A));
L=-tril(A,-1);
U=-A + D -L;
B = D\(L+U);
R = max(abs(eig(B)));
if R <1
n = length(b);
x0 = zeros(n,1);
for i=1:n
    x = x0;
    x(i)=(b(i)-A(i,1:i-1)*x0(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);
end
if norm(x - x0) < TOL
   disp(x)
else
    while norm(x - x0) > TOL
        x0 = x;
        for i=1:n
            x(i)=(b(i)-A(i,1:i-1)*x0(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);
        end
    end
end
else
    x = '雅克比迭代不收敛';
end
end

高斯-赛德尔（G-S）迭代相关的MATLAB函数

function [x] = GS(A,b,TOL)
D=diag(diag(A));
L=-tril(A,-1);
U=-A + D -L;
B = (D-L)\U;
R = max(abs(eig(B)));
if R < 1
n = length(b);
x0 = zeros(n,1);
for i=1:n
    x = x0;
    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
if norm(x - x0) < TOL
   disp(x)
else
    while norm(x - x0) > TOL
        x0 = x;
        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
    end
end
else
    x = '高斯-赛德尔迭代不收敛';
end
end

lagrange差值函数相关的MATLAB函数

function y=lagrange(x0,y0,x)
n=length(x0);
m=length(x);
for i=1:m 
    z=x(i); 
    s=0.0; 
    for k=1:n 
        p=1.0; 
        for j=1:n 
            if j~=k 
                p=p*(z-x0(j))/(x0(k)-x0(j)); 
            end
        end
        s=p*y0(k)+s; 
    end
    y(i)=s; 
end
end

牛顿插值法相关的MATLAB函数
请注意：此处C是第一列为x的值，第二列为y的值2*n的矩阵，n取决于x和y的数量。

function y = Newton_Inter(x,n,C)
%x为插值点，n为插值次数，C为节点的X和Y形成的矩阵
% 判断均插点是否足够
% 以下为计算均插表的代码
[r,c] = size(C);
for j= 3:r+1
    for i=j-1:r
        C(i,j)=( C(i,j-1)-C(i-1,j-1) )/( C(i,1)-C(i+2-j,1));
    end
end
disp(C)
% 核心计算的代码
y=C(1,2);
for i=2:n+1
    S=1;
    for j=i-1:-1:1
        S=S*(x-C(j,1));
    end
y=y+C(i,i+1)*S;
end
disp(y);
end

三次hermit插值法相关的MATLAB函数
复杂函数表示方法如下：

function z=f(a,b,x) 
%b为每一点导数值，a第一行为自变量，第二行为函数值 
	n=size(a,2); 
	p=a(1,:); 
	q=a(2,:); 
	for k=1:n-1 
		if a(1,k)<=x&x<=a(1,k+1)
			z=((x-p(k+1))/(p(k)-p(k+1)))^2*(1+2*(x-p(k))/(p(k+1)-p(k)))*q(k)+... ((x-p(k))/(p(k+1)-p(k)))^2*(1+2*(x-p(k+1))/(p(k)-p(k+1)))*q(k+1)+... 
((x-p(k+1))/(p(k)-p(k+1)))^2*(x-p(k))*b(k)+((x-p(k))/(p(k+1)-p(k)))^2*... 
(x-p(k+1))*b(k+1); 
		break 
		end 
	end
end

直接使用方法：

p = pchip(x,y,xq) 

注意：该式返回与 xq 中的查询点对应的插值 p 的向量。p 的值由 x 和 y 的保形分段三次插值确定。

pp = pchip(x,y) 

注意：该式返回一个分段多项式结构体以用于 ppval 和样条实用工具 unmkpp。

三次样方插值法相关的MATLAB函数
复杂函数表述方法如下：

function [D,h,A,g,M]=threesimple(X,Y)
     n=length(X); 
     A=zeros(n,n);A(:,1)=Y';D=zeros(n-2,n-2);g=zeros(n-2,1);
     for  j=2:n
        for i=j:n
            A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
        end
     end

     for i=1:n-1
         h(i)=X(i+1)-X(i);
     end
     for i=1:n-2
         D(i,i)=2;
         g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2));
     end
     for i=2:n-2
         u(i)=h(i)/(h(i)+h(i+1));
         n(i-1)=h(i)/(h(i-1)+h(i));
         D(i-1,i)=n(i-1);
         D(i,i-1)=u(i);             
     end
     M=D\g;
     M=[1;M;1];         
end
%以下函数需要调用上述函数
function s=threesimple1(X,Y,x)
%自然边界条件函数 
%s函数表示三次样条插值函数插值点对应的函数值
%根据三次样条参数函数求出的D,h,A,g,M
%x表示求解插值点函数点，X为已知插值点        
     [D,h,A,g,M]=threesimple(X,Y);
     n=length(X); m=length(x);    
     for t=1:m
        for i=1:n-1
           if (x(t)<=X(i+1))&&(x(t)>=X(i))
              p1=M(i,1)*(X(i+1)-x(t))^3/(6*h(i));
              p2=M(i+1,1)*(x(t)-X(i))^3/(6*h(i));
              p3=(A(i,1)-M(i,1)/6*(h(i))^2)*(X(i+1)-x(t))/h(i);
              p4=(A(i+1,1)-M(i+1,1)/6*(h(i))^2)*(x(t)-X(i))/h(i);
              s(t)=p1+p2+p3+p4; 
              break;
           else
               s(t)=0; 
           end
        end
     end
end

直接使用方法：

s = spline(x,y,xq)

返回与 xq 中的查询点对应的插值 s 的向量。s 的值由 x 和 y 的三次样条插值确定。

pp = spline(x,y) 

返回一个分段多项式结构体以用于 ppval 和样条实用工具 unmkpp。

上述函数是差值法的MATLAB函数