%Alternative牛顿插值法via541
function [A,y]= newtonzi(X,Y,x)
% Newton插值函数
% X为已知数据点的x坐标
% Y为已知数据点的y坐标
% x为插值点的x坐标
% 函数返回A差商表
% y为各插值点函数值
n=length(X); m=length(x);
for t=1:m
z=x(t); A=zeros(n,n);A(:,1)=Y';
s=0.0; y=0.0; c1=1.0;
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
C=A(n,n);
for k=1:n
p=1.0;
for j=1:k-1
p=p*(z-X(j));
end
s=s+A(k,k)*p;
end
ss(t)=s;
end
y=ss;
A=[X',A];
end

%Hermite插值
function f = Hermite( x,y,y_1,x0 )
%Hermite插值函数
% x为已知数据点的x坐标
% y为已知数据点的y坐标
% y_1为数据点y值导数
% x0为插值点的x坐标
syms t;
f = 0.0;
if(length(x) == length(y))
if(length(x) == length(y_1))
n = length(x);
else
disp('y和y的导数维数不相等');
renturn;
end
else
disp('x和y的维数不相等');
return;
end
%以上为输入判断和确定“n”的值
for i = 1:n
h = 1.0;
a = 0.0;
for j = 1:n
if(j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2);%求得值为(li(x))^2
a = a + 1/(x(i)-x(j)); %求得ai（x）表达式之中的累加部分
end
end
f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
if(i == n)
if(nargin == 4)
f = subs(f, 't' , x0); %输出结果
else
f = vpa(f,6); %输出精度为有效数字为6位的函数表达式
end
end
end

%Hermite插值 主程序
x = [1 1.2 1.4 1.6 1.8];
y = [1 1.0954 1.1832 1.2649 1.3416];
y_1 = [0.5000 0.4564 0.4226 0.3953 0.3727];
disp('显示Hermite插值多项式')
f = Hermite(x,y,y_1)
disp('显示Hermite在插值点的值')
f = Hermite(x,y,y_1,1.8)

%拉格朗日插值
function y=lagrange(x0,y0,x)
% 给定一系列点x0,y0
% x是我们要预测的值，由于可以有多个，因此用向量表示
% y返回我们的估计值，由于可以有多个，因此用向量表示
n = length(x);% 要预测的个数
y = zeros(n);% 初始化，并赋初值0
for k = 1:length(x0)
j_no_k=find((1:length(x0))~=k);% 在这里，find函数用于返回一个向量中不为下标k的元素（下标从1开始）
y1=1;
for j = 1:length(j_no_k)
y1 = y1.*(x-x0(j_no_k(j))); % ∏(x-xj)
end
y = y + y1*y0(k)/prod(x0(k)-x0(j_no_k));% prod函数返回数组元素的乘积
end
end

%三次样条插值
function s=threesimple1(X,Y,x,y0,yn)
[D,h,A,g,M]=three2(X,Y,y0,yn);
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
function [D,h,A,g,M]=three2(X,Y,y0,yn)
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;
if (i==1)
g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2))-h(i)/(h(i)+h(i+1))*y0;
elseif (i==(n-2))
g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2))-(1-h(i)/(h(i)+h(i+1)))*yn;
else
g(i,1)=(6/(h(i+1)+h(i)))*(A(i+2,2)-A(i+1,2));
end
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=[y0;M;yn];
end

%三次样条插值 主程序/画图
x=[0 1 2 3];
y=[3 5 4 1];
y0=1;
yn=1;
x0=0:0.1:3;
s=threesimple1(x,y,x0,y0,yn);
s
plot(x0,s)
hold on
grid on
plot(x,y,'o')
xlabel('自变量 X'), ylabel('因变量 Y')
title('三次样条函数')
legend('三次样条插值点坐标','插值点')

%二分法
a=-2;b=1;
eps=1e-3;
x0=(a+b)/2;
while abs(a-b)>=2*eps
c=(a+b)/2;
if f(a)*f(c)<=0
b=c;
else
a=c;
end
x0=(a+b)/2;
end
disp(['二分法求得的解是：x=',num2str(x0)]);
function y=f(x)
y=sin(x)-6*x-5;
end

%不动点迭代法
g = @(x) fun2(x);
% Initialization
x0 = 0;
tol = 1e-5;
maxIter = 40;
% Test biSection function
[xStar,xRoot] = fixPoint(g,x0,tol,maxIter);
fprintf('The fixed point is: %d\n',xStar);
fprintf('The root of the equation is: %d\n',xRoot);
function [xStar,xRoot] = fixPoint(fun2,x0,tol,maxIter)
% Inputs:
% fun2: a function handle, standing for the function written above
% x0: the initial guess of the fixed point
% tol: the tolerance within which the program can stop
% maxIter: the maximum number of iterations the program is allowed to run
% Outputs:
% xStar: the numerical value of the fixed point
% xRoot: the numerical value of the root
x = zeros(maxIter,1);
x(1) = fun2(x0);
i = 1;
while abs(fun2(x(i))-x(i))>tol && i<maxIter
x(i+1) = fun2(x(i));
xStar = x(i);
i = i+1;
end
xRoot = x(i);
end
function gx = fun2(x)
gx = 1+0.5*sin(x);
end

%牛顿法
function y = Newton(fun_a,fun_b,err,x)
x1 = x - fun_a(x)/fun_b(x);
while ( abs(x-x1)>err )
x=x1;
x1 = x - fun_a(x)/fun_b(x);
end
y = x1
sprintf('牛顿法：结果为:%f',x1)
end

%牛顿弦截法
function [result,tol]=xianjie(f,x,n) %f为函数句柄，x为初始弦截点，n为迭代次数
if nargout == 1
flag=1;
elseif nargout == 2
flag=2;
end
x1=x(1);
x2=x(2);
i=1;
while i < n
x3=x2-f(x2)*(x2-x1)/(f(x2)-f(x1));
if abs(f(x3)) < 1e-8
if flag==1
result=x3;
return;
elseif flag==2
result=x3;
tol=abs(f(x3));
return;
end
end
x1=x2;
x2=x3;
i=i+1;
end
end

%高斯消去
function Gauss(A,b,n)
x=zeros(n,1);
for k=1:n-1
if A(k,k)==0
error('Error')
end
for i=k+1:n
Aik=A(i,k)/A(k,k);
for j=k:n
A(i,j)=A(i,j)-Aik*A(k,j);
end
b(i)=b(i)-Aik*b(k);
end
end
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
s=b(k);
for j=k+1:n
s=s-A(k,j)*x(j);
end
x(k)=s/A(k,k);
end
disp("高斯消去法求解为：")
disp(x)
end

%Jacobi 迭代
function Jacobi(A,b,n)
U = triu(A,1);
L = tril(A,1);
D = tril(triu(A,0),0);
p = eig(D\(L+U));
if max(abs(p)) >= 1
disp("Jacobi迭代不收敛，无法计算")
return
end
x = ones(n,1);
t = 0;
while t<=10 %
for i = 1:n
x1 = x;
for j = 1:n
if j ~= i
x(i) = x(i) + A(i,j) * x1(j);
end
end
x(i) = -(x(i) + b(i)) / A(i,i);
end
t = t + 1;
end
disp("Jacobi迭代求解为：")
disp(x)

%G-S迭代
function GS(A,b,n)
x0 = zeros(n,1);
x1 = zeros(n,1);%
e0 = 1e-3;
M = 100;
m=0;
U = triu(A,1);
L = tril(A,1);
D = tril(triu(A,0),0);
p = eig(D\(L+U));
if max(abs(p)) >= 1
disp("G-S迭代不收敛，无法计算")
return
end
while(1)
for i=1:n
temp=0;
for j=1:i-1
temp=temp+A(i,j)*x1(j);
end
for j=i+1:n
temp=temp+A(i,j)*x0(j);
end
x1(i)=-(temp-b(i))/A(i,i);
end
m=m+1;
if(norm(x1-x0)<e0)
disp("G-S迭代求解为：")
disp(x1)
break;
elseif(m>M)
disp("G-S迭代：达到最大循环次数");
disp(x1);
break;
else
x0=x1;
end
end
end