1、精選優質文檔-傾情為你奉上首先說這個程序并不完善，為了實現通用（1，2，n）格式解題，以及為調用追趕法程序，沒有針對節點數在三個以下的情況進行分類討論。希望能有朋友給出更好的方法。首先，通過函數 sanwanj得到方程的系數矩陣，即追趕法方程的四個向量參數，接下來調用追趕法（在intersanwj函數中），得到三次樣條分段函數系數因子，然后進行多項式合并得到分段函數的解析式，程序最后部分通過判斷輸入值的區間自動選擇對應的分段函數并計算改點的值。附：追趕法程序 chase%function newv,w,newu,newd=sanwj(x,y,x0,y0,y1a,y1b)% 三彎矩樣條插值% 將

2、插值點分兩次輸入，x0 y0 單獨輸入% 邊值條件a的二階導數 y1a 和b的二階導數 y1b，這里建議將y1a和y1b換成y2a和y2b，以便于和三轉角代碼相區別n=length(x);m=length(y);if m=nerror('x or y 輸入有誤，再來');endv=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-1h(k)=x(k+1)-x(k);endv(1)=h0/(h0+h(1);u(1)=1-v(1);d(1)=6*(y(

3、2)-y(1)/h(1)-(y(1)-y0)/h0)/(h0+h(1);%for k=2:n-1v(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);endnewv=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

4、(x);m=length(y);if m=nerror('x or y 輸入有誤，再來');end%重新定義hh=zeros(n,1);h(1)=x(1)-x0;for k=2:nh(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

5、-(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:nfprintf('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);endS=S.'disp(S);fprintf('以上為樣條函數

6、(三彎矩)解析式，顯示為手寫如下：n');pretty(S);%第二部分%是否繼續運行程序myloop=input('繼續運行程序輸入“1”，否則輸入“0”n');if myloopwhile myloopxi=input('輸入需要計算的點的值，并按回車鍵n');if xi>x0|xi<x(n)fprintf('現在開始計算輸入點的插值函數值n');elsefprintf('輸入數值不在插值范圍內，請重新輸入n');xi=input('輸入需要計算的點的值，并按回車鍵n');end% 確定輸入

7、的數值應該使用哪個解析式newx=x0;x;r,suoy=min(abs(newx-xi);fprintf('輸入點的插值函數值為:nn');fprintf('t');if xi<=newx(suoy)f=subs(S(suoy-1),X,xi);elsef=subs(S(suoy),X,xi);enddisp(f);myloop=input('繼續計算輸入“1”，終止計算輸入“0”n');endelsereturn;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-1v(i)=c(i)/(b(i)-v(i-1)*a(i-1);u(i)=(d(i)-u(i-1)*a(i-1)/(b(i)-v(