sanwj.rar_三次样条拟合_导数_导数法_插值_点拟合

%%%%%%%%%%%%%% function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b) % 三弯矩样条插值 % 将插值点分两次输入，x0 y0 单独输入 % 边值条件a的二阶导数 y1a 和b的二阶导数 y1b，这里建议将y1a和y1b换成y2a和y2b，以便于和三转角代码相区别 n=length(x);m=length(y); if m~=n error('x or y 输入有误，再来'); end v=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1); w=2*ones(n+1); h0=x(1)-x0; h=zeros(n-1,1); for k=1:n-1 h(k)=x(k+1)-x(k); end v(1)=h0/(h0+h(1)); u(1)=1-v(1); d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1)); % for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k)); u(k)=1-v(k); d(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k)); end newv=[v;1]; newu=[1;u]; d0=6*((y(1)-y0)/h0-y1a)/h0; d(n)=6*(y1b-(y(n)-y(n-1))/h(n-1))/h(n-1); newd=[d0;d]; %%%%%%%%%%%% function intersanwj(x,y,x0,y0,y1a,y1b) % 三弯矩样条插值 %第一部分 n=length(x);m=length(y); if m~=n error('x or y 输入有误，再来'); end %重新定义h h=zeros(n,1); h(1)=x(1)-x0; for k=2:n h(k)=x(k)-x(k-1); end %sptep1 调用三弯矩函数 [a,b,c,d]=sanwj(x,y,x0,y0,y1a,y1b); % 三对角方程 M=chase(a,b,c,d); % 求插值函数 fprintf('三次样条(三弯矩)插值的函数表达式\n'); syms X ; fprintf('S0--1:\n'); S(1)=collect(((1/6)*M(2)*(X-x0).^3-(1/6)*M(1)*(X-x(1)).^3+(y(1)-(M(2)*h(1).^2)/6)*(X-x0)-(y0-(M(1)*h(1).^2)/6)*(X-x(1)))/h(1)); for k=2:n fprintf('S%d--%d:\n',k-1,k); S(k)=collect(((1/6)*M(k+1)*(X-x(k-1)).^3-(1/6)*M(k)*(X-x(k)).^3+(y(k)-(M(k+1)*h(k).^2)/6)*(X-x(k-1))-(y(k-1)-(M(k)*h(k).^2)/6)*(X-x(k)))/h(k)); end S=S.'; disp(S); fprintf('以上为样条函数(三弯矩)解析式，显示为手写如下：\n'); pretty(S); %第二部分 %是否继续运行程序 myloop=input('继续运行程序输入“1”，否则输入“0”\n'); if myloop while myloop xi=input('输入需要计算的点的值，并按回车键\n'); if xi>x0|xi<x(n) fprintf('现在开始计算输入点的插值函数值……\n'); else fprintf('输入数值不在插值范围内，请重新输入\n'); xi=input('输入需要计算的点的值，并按回车键……\n'); end % 确定输入的数值应该使用哪个解析式 newx=[x0;x]; [r,suoy]=min(abs(newx-xi)); fprintf('输入点的插值函数值为:\n\n');fprintf('\t'); if xi<=newx(suoy) f=subs(S(suoy-1),X,xi); else f=subs(S(suoy),X,xi); end disp(f); myloop=input('继续计算输入“1”，终止计算输入“0”\n'); end else return; end %%%%%%%%%%%% function [x]=chase(a,b,c,d) %追赶法解性方程组 a是下三角b是对角线c是上三角 d是常数项 %输入的a b c d 均为列向量 n=length(b); u=zeros(n,1); v=zeros(n,1); x=zeros(n,1); %追 v(1)=c(1)/b(1);u(1)=d(1)/b(1); for i=2:n-1 v(i)=c(i)/(b(i)-v(i-1)*a(i-1)); u(i)=(d(i)-u(i-1)*a(i-1))/(b(i)-v(i-1)*a(i-1)); end u(n)=(d(n)-u(n-1)*a(n-1))/(b(n)-v(n-1)*a(n-1)); %赶 x(n)=u(n); for i=n-1:-1:1 x(i)=u(i)-v(i)*x(i+1); end