数值分析程序代码(MATLAB)_牛顿插值法_三次样条插值多项式_插值型求积_多项式插值_样条插值_

function [p,M,S]=Spline3(flag) % Made by 付林 % Input % flag 边界条件类型，flag==1(固支边界)，flag==2(简支边界) % f为输入函数 f=1/(1+x^2) x=[-5,5] % Output % p 输出多项式 % X为给定的差值节点 X=linspace(-5,5,11); syms x; % syms x; f=1/(1+x^2) f=1/(1+x^2); % Y为给定差值节点X对应的函数值 Y=1./(1+X.^2); lX=length(X); % 构造临时数组r,u,d,h r=[];u=[];d=[];h=[]; for i=1:lX-1 % 求解每个区间长度h(i) h(i)=X(i+1)-X(i); end %% 固支边界条件（I型边界条件) if (flag==1) % 判断取的是什么边界条件，固支边界条件 %% 步骤1、2 求解参数r(i),u(i),d(i)过程，并补充两边界条件 for j=2:lX-1 % 求解r(i),为了三弯矩求解方便，故这里设置为u(i) u(j-1)=h(j)/(h(j-1)+h(j)); % 求解u(i),为了三弯矩求解方便，故这里设置为r(i) r(j-1)=h(j-1)/(h(j-1)+h(j)); % 求解d(i) d(j-1)=3*(u(j-1)*(Y(j)-Y(j-1))/h(j-1)+r(j-1)*(Y(j+1)-Y(j))/h(j)); end % 对函数f求一阶导,并构造为内联函数df1 df1=inline(diff(f),'x'); % S'(X0)=f'(X0) m(1)=feval(df1,X(1)); % S'(Xn)=f'(Xn) m(2)=feval(df1,X(lX)); p=[r;u;d]; % 将三弯矩方程化为追赶法标准方程 d(1)=d(1)-u(1)*m(1); % 将三弯矩方程化为追赶法标准方程 d(lX-2)=d(lX-2)-r(lX-2)*m(2); %% 步骤3 追赶法求解m参数 for i=1:length(d) % 对角三角矩阵的对角线 b(i)=2; end % 对角线下的数组 a=u;a(1)=0; % 对角线上的数组 c=r;c(lX-2)=0; % 调用ZhuiGanFa()函数，求解M的解(但不包括M0和Mn) M=ZhuiGanFa(a,b,c,d); M=[m(1) M m(2)]; %% 步骤4 回代到方程组S(x)中 for i=1:length(h) S1=Y(i)*((x-X(i+1))^2*(h(i)+2*(x-X(i))))/(h(i))^3; S2=Y(i+1)*((x-X(i))^2*(h(i)-2*(x-X(i+1))))/(h(i))^3; S3=M(i)*((x-X(i))*(x-X(i+1))^2)/(h(i))^2; S4=M(i+1)*((x-X(i+1))*(x-X(i))^2)/(h(i))^2; % 简化多项式S(x) S(i,1)=simplify(S1+S2+S3+S4); end %% 简支边界条件（II型边界条件） elseif (flag==2) % 判断取的是什么边界条件，简支边界条件 %% 步骤1、2 求解参数r(i),u(i),d(i)过程，并补充两边界条件 for j=2:lX-1 % 求解r(i) r(j-1)=h(j)/(h(j-1)+h(j)); % 求解u(i) u(j-1)=h(j-1)/(h(j-1)+h(j)); % 求解d(i) d(j-1)=(6/(h(j-1)+h(j)))*((Y(j+1)-Y(j))/h(j)+(Y(j)-Y(j-1))/h(j-1)); end % 对函数f求二阶导,并构造为内联函数df2 df2=inline(diff(f,2),'x'); % S''(X0)=f''(X0) m(1)=feval(df2,X(1)); % S''(Xn)=f''(Xn) m(2)=feval(df2,X(lX)); p=[r;u;d]; % 将三弯矩方程化为追赶法标准方程 d(1)=d(1)-u(1)*m(1); % 将三弯矩方程化为追赶法标准方程 d(lX-2)=d(lX-2)-r(lX-2)*m(2); %% 步骤3 追赶法求解m参数 for i=1:length(d) % 对角三角矩阵的对角线 b(i)=2; end % 对角线下的数组 a=u;a(1)=0; % 对角线上的数组 c=r;c(lX-2)=0; % 调用ZhuiGanFa()函数，求解M的解(但不包括M0和Mn) M=ZhuiGanFa(a,b,c,d); M=[m(1) M m(2)]; %% 步骤4 回代到方程组S(x)中 for i=1:length(h) S1=M(i)*(X(i+1)-x)^3/(6*h(i)); S2=M(i+1)*(x-X(i))^3/(6*h(i)); S3=(Y(i)-(M(i)*(h(i))^2)/6)*(X(i+1)-x)/h(i); S4=(Y(i+1)-(M(i+1)*(h(i))^2)/6)*(x-X(i))/h(i); % 简化多项式S(x) S(i,1)=simplify(S1+S2+S3+S4); end end %% 绘制图像三次样条插值图像 for k=1:length(S) fh=inline(S(k,1),'x'); fplot(fh,[X(k),X(k)+h(k)]); hold on; end %% 绘制原始图像 f=inline(f,'x'); fplot(f,[-5,5],'k'); % 绘制图例 legend('三次样条','原始图像');