【導讀】,要求誤差小于610?,給出實現(xiàn)算法。具體的要求參見所附的相關文檔。針對本專業(yè)中所碰到的實際問題，提煉一。個使用方程組進行求解的例子，并對求解過程進行分析、求解。以及2020年的產(chǎn)量?；囊?guī)律,試計算這一天的平均氣溫。

　　

【正文】 程組函數(shù) %求解三對角線性方程組系數(shù)存放在 a， b， c 中， d 為右端向量 function [x] = TSS_wlg(a,b,c,d,n) u=zeros(1,n)。 y=zeros(1,n)。 l=zeros(1,n1)。 u(1)=b(1)。 y(1)=d(1)。 for k=2:n l(k1)=a(k1)/u(k1)。 u(k)=b(k)l(k1)*c(k1)。 y(k)=d(k)l(k1)*y(k1)。 end x(n)=y(n)/u(n)。 for k=n1:1:1 x(k)=(y(k)c(k)*x(k+1))/u(k)。 end 三次樣條函數(shù) %三次樣條方法 % flag 表端點邊界條件類型： % flag=0 : 自然樣條 (端點二階導數(shù)為 0)； % flag=1 : 第一類邊界條件 (端點一階導數(shù)給定 )； % flag=2 : 第二類邊界條件 (端點二階導數(shù)給定 )； % y00,ynn 表左右端點處的在邊界條件值。 function [M] =yangjiao_wlg( x,y,flag,y00,ynn ) n=length(x)1。 M=zeros(1,n+1)。 h=zeros(1,n)。 a=zeros(1,n)。 b=zeros(1,n+1)。 c=zeros(1,n)。 計算方法實驗報告 37 for i=1:n+1 M(i)=y(i)。 end for k=1:2 for i=n+1:1:k+1 M(i)=(M(i)M(i1))/(x(i)x(ik))。 end end h(1)=x(2)x(1)。 for i=2:n h(i)=x(i+1)x(i)。 c(i)=h(i)/(h(i)+h(i1))。 a(i1)=1c(i)。 b(i)=2。 M(i)=6*M(i+1)。 end b(1)=2。b(n+1)=2。 % 自然樣條端點條件 (端點二階導數(shù)為零 ) % if flag==0 M(1)=0。M(n+1)=0。c(1)=0。a(n)=0。 end % % % 端點一階導數(shù)條件 % if flag==1 d(1)=6/(x(2)x(1))*((y(2)y(1))/(x(2)x(1))y00)。 d(n)=6/(x(n)x(n1))*(ynn(y(n)y(n1))/(x(n)x(n1)))。 c(1)=1。 a(n)=1。 end % % % 端點二階導數(shù)條件 % if flag==2 d(1)=12*(x(2)x(1))*.... ((((y(4)y(3))/(x(4)x(3))(y(3)y(2))/(x(3)x(2)))/(x(4)x(2))).... (((y(3)y(2))/(x(3)x(2))(y(2)y(1))/(x(2)x(1)))/(x(3)x(1))))/(x(4)x(1))。 d(n)=12*(x(n)x(n1))*.... ((((y(n)y(n1))/(x(n)x(n1))(y(n1)y(n2))/(x(n1)x(n2)))/(x(n)x(n2))).... (((y(n1)y(n2))/(x(n1)x(n2))(y(n2)y(n3))/(x(n2)x(n3)))/(x(n1)x(n3))))/(x(n)x(n3))。 c(1)=2。 a(n)=2。 end % % 計算方法實驗報告 38 M=TSS_wlg(a,b,c,M,n+1)。 查找插值點函數(shù) %查找給定點 x0 在的位置，以 x(k)下表 k 作為查詢結(jié)果 function k= findk_wlg( x,y,xi) n=length(x)。 k=1。 for i=2:n if xix(i) k=i。 break。 else k=i+1。 end end 計算方法實驗報告 39 附錄五 問題 5Matlab 求解 法方程解法程序 %法方程法最小二乘法 clc clear % x=input(39。請輸入橫坐標 x=:39。) % y=input(39。請輸入縱坐標 y=:39。) x=0:24。 y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31.... 29 27 25 24 22 20 18 17 16]。 %依次使用二次函數(shù)、三次函數(shù)、四次函數(shù)的最小二乘法擬合 %利用法方程法，其中 A=G39。*G， b=G39。*y for m=2:4 A=zeros(m+1,m+1)。 for i=0:m for j=0:m A(i+1,j+1)=sum(x.^(i+j))。 end bb(i+1)=sum(x.^i.*y)。 end aa=A\bb39。 %左右翻轉(zhuǎn) disp(39。擬合多項式如下及相應的系數(shù) 39。) p=fliplr(aa39。) Px=poly2sym(p) xi=0::24。 yi=polyval(p,xi)。 y1=polyval(p,x)。 disp(39。估算誤差 39。) e=sum((yy1).^2) subplot(2,2,m1)。 plot(x,y,39。.39。,xi,yi,39。r39。)。 title(39。duoxiangshi39。)。 end %指數(shù)型的函數(shù) C=a*exp(b*(tc)^2)最小二乘法擬合 x2=x。 y2=log(y)。 m=2。 B=zeros(m+1,m+1)。 計算方法實驗報告 40 for i=0:m for j=0:m B(i+1,j+1)=sum(x2.^(i+j))。 end cc(i+1)=sum(x2.^i.*y2)。 end z=B\cc39。 p=fliplr(z39。) b=p(1) c=p(2)/(2*b) a=exp(p(3)+b*c^2) xi=0::24。 yi=exp(polyval(p,xi))。 y1=exp(polyval(p,x2))。 disp(39。估算誤差 39。) e=sum((yy1).^2) subplot(2,2,4)。 plot(x2,y,39。.39。,xi,yi,39。r39。)。 title(39。zhishu39。)。 正交分解法主程序 %正交方法最小二乘法 clc clear %x=input(39。請輸入橫坐標 x=:39。) %y=input(39。請輸入縱坐標 y=:39。) x=0:24。 y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31.... 29 27 25 24 22 20 18 17 16]。 %依次使用二次函數(shù)、三次函數(shù)、四次函數(shù)的最小二乘法擬合 %二次函數(shù) % disp(39。擬合二次函數(shù)如下及相應的系數(shù) 39。) z=G_wlg(x,y,2)。 zz=QG_wlg(z)。 %[q,r]=qr(z) aa=jief_wlg(zz)。 p1=fliplr(aa39。) Px1=poly2sym(p1) xi=0::24。 yi=polyval(p1,xi)。 y1=polyval(p1,x)。 disp(39。估算誤差 39。) 計算方法實驗報告 41 e1=wucha_wlg(zz) subplot(2,2,1)。 plot(x,y,39。.39。,xi,yi,39。r39。)。 title(39。二次函數(shù) 39。)。 % %三次函數(shù) % disp(39。擬合三次函數(shù)如下及相應的系數(shù) 39。) z=G_wlg(x,y,3)。 zz=QG_wlg(z)。 %[q,r]=qr(z) aa=jief_wlg(zz)。 p2=fliplr(aa39。) Px2=poly2sym(p2) xi=0::24。 yi=polyval(p2,xi)。 y1=polyval(p2,x)。 disp(39。估算誤差 39。) e2=wucha_wlg(zz) subplot(2,2,2)。 plot(x,y,39。.39。,xi,yi,39。r39。)。 title(39。三次函數(shù) 39。)。 % %四次函數(shù) % disp(39。擬合四次函數(shù)如下及相應的系數(shù) 39。) z=G_wlg(x,y,4)。 zz=QG_wlg(z)。 %[q,r]=qr(z) aa=jief_wlg(zz)。 p3=fliplr(aa39。) Px3=poly2sym(p3) xi=0::24。 yi=polyval(p3,xi)。 y1=polyval(p3,x)。 disp(39。估算誤差 39。) e3=wucha_wlg(zz) subplot(2,2,3)。 plot(x,y,39。.39。,xi,yi,39。r39。)。 title(39。四次函數(shù) 39。)。 % %指數(shù)型的函數(shù) C=a*exp(b*(tc)^2)最小二乘法擬合 % disp(39。擬合指數(shù)函數(shù)的系數(shù)如下 39。) 計算方法實驗報告 42 x2=x。 y2=log(y)。 z=G_wlg(x2,y2,2)。 zz=QG_wlg(z)。 %[q,r]=qr(z) aa=jief_wlg(zz)。 p4=fliplr(aa39。) b=p4(1) c=p4(2)/(2*b) a=exp(p4(3)+b*c^2) xi=0::24。 yi=exp(polyval(p4,xi))。 y1=exp(polyval(p4,x2))。 disp(39。估算誤差 39。) e4=wucha_wlg(zz) subplot(2,2,4)。 plot(x2,y,39。.39。,xi,yi,39。r39。)。 title(39。指數(shù)函數(shù) 39。)。 % G 矩陣求解函數(shù) %G 矩陣 function [G] = G_wlg(x,y,n) %得到 G 矩陣，最后一列存放 y 值， n 表示擬合次數(shù)， 1,x,x^2... m=length(x)。 G=zeros(m,n+2)。 for i=1:m for j=1:n+1 G(i,j)=(x(i))^(j1)。 end end for i=1:m G(i,n+2)=y(i)。 end RO 與 h h2 矩陣函數(shù) %RO 與 h h2 矩陣 function [GG] = QG_wlg(G) %依次變換得到 Gk， Qk，最終得到 RO 與 h h2 [m,n]=size(G)。 n=n1。 w=zeros(m,1)。 計算方法實驗報告 43 for k=1:n %依次變換得到 Qk gik=0。 for i=k:m gik=gik+G(i,k)^2。 end dlt=(sign(G(k,k)))*sqrt(gik)。 w(k)=G(k,k)dlt。 for j=k+1:m w(j)=G(j,k)。 end blt=dlt*w(k)。 %依次變換得到 Gk G(k,k)=dlt。 %k 取值 for j=k+1:n+1 wg=0。 for ii=k:m wg=wg+w(ii)*G(ii,j)。