【正文】 clear disp(39。) A=input(39。) %判斷被分解矩陣的維數(shù) [n,n]=size(A)。 L(j,1)=A(j,1)/U(1,1)。 for k=2:n for j=k:n for m=1:(k1) %計算公式中對應(yīng)的 l*u 的和 uu(m)=L(k,m)*U(m,j)。 for m=1:(k1) uu(m)=L(j,m)*U(m,k)。 end end L U 計算方法實驗報告 30 附錄三 問題 3Matlab 求解 主程序： %列主元高斯消去法求解方程解 clc。 %打開二進制文件 [fid message]=fopen(39。)。 else %id:為該數(shù)據(jù)文件的標(biāo)識，值為 F1E1D1A0 id=fread(fid,1,39。)。int3239。 if (ver==257)||(ver==258) ver=dec2hex(ver) %n:表示方程的階數(shù) n=fread(fid,1,39。) %系數(shù)矩陣 ab=fread(fid,[n n],39。)。 whos ab %右端系數(shù) d=fread(fid,[n 1],39。) whos d %列主元高斯消去法 gauss(Ab,d) else if ver==513 ver=dec2hex(ver) %n:表示方程的階數(shù) n=fread(fid,1,39。) %系數(shù)矩陣 %為下帶寬 p=fread(fid,1,39。) %為上帶寬 計算方法實驗報告 31 q=fread(fid,1,39。) %系數(shù)矩陣 ab=fread(fid,[p+q+1,n],39。)。 % %驗證，轉(zhuǎn)換為一般系數(shù)矩陣 %abb=zeros(n,n)。 %end %end %abb %驗證，轉(zhuǎn)換為一般系數(shù)矩陣 % %列主元高斯消去法 %gauss(abb,d) whos abbb %右端系數(shù) d=fread(fid,[n 1],39。)。 求解一般系數(shù)矩陣的列主元高斯消去法函數(shù)： %列主元高斯消去法 function gauss(A,b) [m,n]=size(A)。 for k=1:n1 max=k。 end end %交換主元所對應(yīng)的行 temp=Ab(k,k:n+1)。 Ab(max,k:n+1)=temp。 Ab(i,k+1:n+1)=Ab(i,k+1:n+1)+Ab(i,k)*Ab(k,k+1:n+1)。 x(n,1)=Ab(n,n+1)/Ab(n,n)。 for j=i+1:n sum=sum+x(j,1)*Ab(i,j)。 end %顯示方程組的解 disp(39。) for xxp=1:n fprintf(39。,xxp,x(xxp)) end 求解壓縮系數(shù)矩陣的列主元高斯消去法函數(shù)： %列主元高斯消去法 ,對于壓縮矩陣方程組 %A為壓縮系數(shù)矩陣 %b 為右端矩陣 %p 為下帶寬 %q 為上帶寬 function gausskz(A,b,p,q) [n,m]=size(A)。%p 有關(guān) %增廣陣 Ab=[A c b]。 Ab(k,2:6)=Ab(k+1,1:5)。 %交換右端系數(shù) temp1=Ab(k,7)。 Ab(k+1,7)=temp1。 Ab(k+1,1:5)=Ab(k+1,1:5)+a*Ab(k,2:6)。 end for i=1:n for j=1:p+q+1 Abb(i,j)=Ab(i,j+1)。 end %回代解方程 x=zeros(n,1)。 for i=n1:1:n4 sum=0。 end x(i,1)=(Abb(i,6)sum)/Abb(i,1)。 for jj=ii+1:ii+4 sum1=sum1+x(jj,1)*Abb(ii,jjii+1)。 end %顯示方程組的解 disp(39。) for xxp=1:n fprintf(39。,xxp,x(xxp)) end 計算方法實驗報告 34 附錄四 問題 4Matlab 求解 主程序 %三次樣條插值及 Newton 插值 clc。 x=[1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2020 2020]。 n=length(x)。o39。 %三次樣條插值 yi=zeros(size(xi))。 end yi(111)=y(n)。 for j=112:121 h=x(12)x(11)。 x2=xi(j)x(11)。 yi(j)=(yi(j)+(y(11)M(11)*h*h/6)*x1+(y(12)M(12)*h*h/6)*x2)/h。 disp(39。) [Nx,N0]=Newton_wlg(x,y,2020) disp(39。) [Nx,N01]=Newton_wlg(x,y,2020) [Nx,Ni]=Newton_wlg(x,y,xi1)。 plot(xi1,Ni,39。)。Newton39。 subplot(2,1,2)。.r39。 title(39。)。三次樣條插值 2020 年的產(chǎn)量及插值函數(shù) 39。三次樣條插值 2020 年的產(chǎn)量 39。 %多項式自變量 x syms x。 y1=0。 Ni=1。 end yi(i)=y1(i+1)。 Nx=Nx+Ni*yi(i)。 end %x0 時值 N0=subs(Nx,39。,x0)。 M=zeros(1,n+1)。 k=findk_wlg(x,y,xi)。 x1=x(k)xi。 yi=M(k1)*x1*x1*x1/6+M(k)*x2*x2*x2/6。 求對應(yīng)點值函數(shù)值及該點對應(yīng)插值函數(shù) %求對應(yīng)點的值 function yi=yihanshu_wlg( x,y,xi) syms t。 M=zeros(1,n+1)。 k=findk_wlg(x,y,xi)。 x1=x(k)t。 y0=M(k1)*x1*x1*x1/6+M(k)*x2*x2*x2/6。 f=simplify(y0) yi=subs(y0,39。,xi)。 y=zeros(1,n)。 u(1)=b(1)。 for k=2:n l(k1)=a(k1)/u(k1)。 y(k)=d(k)l(k1)*y(k1)。 for k=n1:1:1 x(k)=(y(k)c(k)*x(k+1))/u(k)。 function [M] =yangjiao_wlg( x,y,flag,y00,ynn ) n=length(x)1。 h=zeros(1,n)。 b=zeros(1,n+1)。 計算方法實驗報告 37 for i=1:n+1 M(i)=y(i)。 end end h(1)=x(2)x(1)。 c(i)=h(i)/(h(i)+h(i1))。 b(i)=2。 end b(1)=2。 % 自然樣條端點條件 (端點二階導(dǎo)數(shù)為零 ) % if flag==0 M(1)=0。c(1)=0。 end % % % 端點一階導(dǎo)數(shù)條件 % if flag==1 d(1)=6/(x(2)x(1))*((y(2)y(1))/(x(2)x(1))y00)。 c(1)=1。 end % % % 端點二階導(dǎo)數(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))。 c(1)=2。 end % % 計算方法實驗報告 38 M=TSS_wlg(a,b,c,M,n+1)。 k=1。 break。 end end 計算方法實驗報告 39 附錄五 問題 5Matlab 求解 法方程解法程序 %法方程法最小二乘法 clc clear % x=input(39。) % y=input(39。) x=0:24。 %依次使用二次函數(shù)、三次函數(shù)、四次函數(shù)的最小二乘法擬合 %利用法方程法，其中 A=G39。*y for m=2:4 A=zeros(m+1,m+1)。 end bb(i+1)=sum(x.^i.*y)。 %左右翻轉(zhuǎn) disp(39。) p=fliplr(aa39。 yi=polyval(p,xi)。 disp(39。) e=sum((yy1).^2) subplot(2,2,m1)。.39。r39。 title(39。)。 y2=log(y)。 B=zeros(m+1,m+1)。 end cc(i+1)=sum(x2.^i.*y2)。 p=fliplr(z39。 yi=exp(polyval(p,xi))。 disp(39。) e=sum((yy1).^2) subplot(2,2,4)。.39。r39。 title(39。)。請輸入橫坐標(biāo) x=:39。請輸入縱坐標(biāo) y=:39。 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ù)如下及相應(yīng)的系數(shù) 39。 zz=QG_wlg(z)。 p1=fliplr(aa39。 yi=polyval(p1,xi)。 disp(39。) 計算方法實驗報告 41 e1=wucha_wlg(zz) subplot(2,2,1)。.39。r39。 title(39。)。擬合三次函數(shù)如下及相應(yīng)的系數(shù) 39。 zz=QG_wlg(z)。 p2=fliplr(aa39。 yi=polyval(p2,xi)。 disp(39。) e2=wucha_wlg(zz) subplot(2,2,2)。.39。r39。 title(39。)。擬合四次函數(shù)如下及相應(yīng)的系數(shù) 39。 zz=QG_wlg(z)。 p3=fliplr(aa39。 yi=polyval(p3,xi)。 disp(39。) e3=wucha_wlg(zz) subplot(2,2,3)。.39。r39。 title(39。)。擬合指數(shù)函數(shù)的系數(shù)如下 39。 y2=log(y)。 zz=QG_wlg(z)。 p4=fliplr(aa39。 yi=exp(polyval(p4,xi))。 disp(39。) e4=wucha_wlg(zz) subplot(2,2,4)。.39。r39。 title(39。)。 G=zeros(m,n+2)。 end end for i=1:m G(i,n+2)=y(i)。 n=n1。 計算方法實驗報告 43 for k=1:n %依次變換得到 Qk gik=0。 end dlt=(sign(G(k,k)))*sqrt(gik)。 for j=k+1:m w(j)=G(j,k)。 %依次變換得到 Gk G(k,k)=dlt。 for ii=k:m wg=wg+w(ii)*G(ii,j