計算幾年前的只做參考切照抄上機報告

本次上機，編程語言采用語言，編譯環(huán)境為7.6.0，運行平臺為WindowsXP。GaussAx=b6464117513273129227111-1Gauss1-2Gauss%列主元gaussfunction %s消去可以進行否，x %矩陣行和列數(shù)m，nifm~=n|rank(a)~=n %a0,消去不能進行，輸出錯誤信息 fork=1:n-1 %m

fork=n-1:-1:1 a=[-64--11-1-275-6-1327--4312129----72-27-b=[2958 %調用Gaussif~s disp('Error!Thereissomethingwrongwiththematrix!'); disp('Thesolutionis:Thesolution71的列向量。將矩陣a與得到的解向量相乘，剛矩陣LDLT分解與Cholesky

Aijij

i,i。 min(i,j)-2

i然后將矩陣的Choleskyfunctions,g]=cholesky(a)，矩陣的LDLT%矩陣的LDLT分解functions,l,d]=ldlt(a)ifa~=a' %b存放矩陣a的對角元素，矩陣D %矩陣anfork=1:nif~b(k) %如果矩陣D0，出現(xiàn)錯誤，停止計算 ifs

矩陣的Cholesky%矩陣的cholesky分解functions,g]=cholesky(a)ifa~=a'|min(eig(a))<=0 %矩陣a維數(shù)nforforj=1:i-a(i,j)=(a(i,j)-a(i,1:j-1)*a(j,1:j-1)')/a(j,j對角線下方元素計算);%

%a%ldlt %ldltif~s disp('Error!Theldlt positioncannotgo!');

disp('TheMatrixLis:');disp(l); %輸出矩陣Ldisp('TheMatrixDis:');disp(d); %輸出矩陣Ddisp('TheMatrixLTis:');disp(l');%輸出矩陣%進行cholesky %調用choleskyif~s disp('Error!Thecholesky positioncannotgo!');

disp('TheMatrixGis:');disp(g); %輸出矩陣Gdisp('TheMatrixGTis: %輸出矩陣20x20的單位矩陣；LT20x201，其余元素為-1Cholesky分解：G20x201，其余元素為-1GT20x201，其余元素為-1的上三角矩陣。通過計算L×D×LT的結果等于a，可知LDLT分解正確；G×GT的結果等于a，可知Cholesky矩陣LDLT分解與Cholesky

20階三對角方程

Tx

的解，其中

x(x,

T

,2,f(3,2, 4

)Tfunctionx=chase(a,b,c,d)是三對角矩陣對角線下方的元素，b為矩陣對角線上的元素，為對角線上方的元素，x返回方程組的解向量。過程分為追和趕兩個過程，具體算法見P49functionx=chase(a,b,c,d)fork=2:n fork=n-1:- d=[3,2+zeros(1,n-2),3];%右端向量disp('Thesolution 201的列向量。Thesolutionis:Jacobi

10i225i3 45i4 5i230i5Jacobi首先將Jacobifunction[s,i,y]=jacobi(a,b,r)，其中s返Gauss-seidel迭代也編寫成子函數(shù)的形式function[s,i,y]=體算法是先算出迭代矩陣G，然后逐步迭代，直至解達到指定精度。主函數(shù)中將矩陣a，向量b和誤差r輸入后，分別調用函數(shù)jacobi(a,b,r)gaussseidel(a,b,r)Jacobi%Jacobi迭代法解方程組functions,i,y]=jacobi(a,b,r) %矩陣aifm~=n|rank(a)~=n

%da%若矩陣a00 %0if %0for fori=1:10^5 ifmax(abs(y-

Gauss-Seidel%Gauss-Seidel迭代法解方程組functions,i,y]=gaussseidel(a,b,r) %矩陣aifm~=n|rank(a)~=n

%da%若矩陣a00 %0if %0for forj=1:n-1 fori=1:10^5 ifmax(abs(y-

- -10-0-10 -00-000000%jacobiif~sdisp('Error!Jacobiitingmethodcannotgo!');elseifi==10^5disp('Error!ItisdivergentwhenusingJacobiitingmethod!'); disp(['After',num2str(i),'timesJacobiition,thesolutionconvergeon:']);%采用Gauss-Seidelif~sdisp('Error!Gauss-seidelitingmethodcannotgo!');elseifi==10^5disp('Error!ItisdivergentwhenusingGauss-seidelitingmethod!'); disp(['After',num2str(i),'timesGauss-seidelition,thesolutionconvergeon:']);After8timesJacobiition,thesolutionconvergeon:After4timesGauss-seidelition,thesolutionconvergeon:由結果可以看出，利用Jacobi8次迭代就可達到要求的精度，而Seidel迭代法收斂較快。

已知f(x)

125

(1x1),

n5,10,20作

f(x)

xi1

2i,(i0,1,2 ,,n),n)

f(xi)

xi

yi

(i0,1,2

插值多項式Nn

Sn(x)3、對n5,20

xk1

yk

f(xk

Nn(xk)

Sn(xk (k1, ,99)4、計算:E(Nnk

Nn(xk

，E(Sn)k

ykSn(xk functions=myspline(m,y)，其中m,y為已知的插值數(shù)據點，s返回插值后的各個區(qū)間的多項式。樣條函數(shù)采用的Mfunctionx=chase(a,b,c,d)a、b、c、d向量中的個別值就可以，具體見P109。new=newton(m,y)和functions=myspline(m,y)進行插值，計算完后將各個插值函Newton%Newton插值多項式functionnew=newton(m,y)fori=2:n functions=myspline(m,y)symsx; fori=1:n- %指定nsymsx %符號數(shù)x %定義函數(shù)f(x) disp(['Whenn=',num2str(n),',thevaluesoff(x)atthepointsofx(i)are:']);%Newton disp(['Whenn=',num2str(n),',N(x)=']);disp(nout);fori=1:n a=sym2poly(s(i));b=poly2sym(a);b=vpa(b,4整理多項式并輸出disp(['Intherangeof',num2str(m(i)),',',num2str(m(i+1)),'],S(x)=']); %計算xk的值 %計算f(xk),N(xk)fori=1:length(xk) % %計算E(Nn),E(Sn),并輸出title('Interpolation');legend('Originalfunction','NewtonInterpolation','Spline') Whenn=5,thevaluesoff(x)atthepointsofx(i) Intherangeof[-1,-0.6],S(x)=1.721*x^3+5.162*x^2+5.040*x+1.638Intherangeof[-0.6,-0.2],S(x)=-3.315*x^3-3.902*x^2-.3978*x+.5500Intherangeof[-0.2,0.2],S(x)=.5765-1.913*x^2Intherangeof[0.2,0.6],S(x)=3.315*x^3-Intherangeof[0.6,1],S(x)=1.638-5.040*x+5.162*x^2-1.721*x^3n=10Whenn=10,thevaluesoff(x)atthepointsofx(i) Intherangeof[-1,-0.8],S(x)=.3417*x^3+1.025*x^2+1.113*x+.4683Intherangeof[-0.8,-0.6],S(x)=.8933*x^3+2.349*x^2+2.172*x+.7507Intherangeof[-0.6,-0.4],S(x)=.8364*x^3+2.246*x^2+2.111*x+.7384Intherangeof[-0.4,-0.2],S(x)=13.41*x^3+17.33*x^2+8.145*x+1.543Intherangeof[-0.2,0],S(x)=-54.47*x^3-23.39*x^2-.2872e-14*x+1.Intherangeof[0,0.2],S(x)=54.47*x^3-23.39*x^2+.2872e-14*x+1.Intherangeof[0.2,0.4],S(x)=-13.41*x^3+17.33*x^2-8.145*x+1.543Intherangeof[0.4,0.6],S(x)=-.8364*x^3+2.246*x^2-2.111*x+.7384Intherangeof[0.6,0.8],S(x)=-.8933*x^3+2.349*x^2-2.172*x+.7507Intherangeof[0.8,1],S(x)=-.3417*x^3+1.025*x^2-1.113*x+.4683中的n修改即可，此處不再贅述。為了比較n取不同值時插值函數(shù)的變化，將n從上圖中可以看出，當n=5時，插值函數(shù)與原函數(shù)有較大的誤差，而當n取Newton從上圖中看出隨著n的增大，插值精度反而降低，特別是在區(qū)間兩端，插值式，而要取相對較小的n用n4

4x35x22x 08(x (10

dt

0,1

y(x)

(x

將積分號里面的函數(shù)離散化，其自變量為t，并用復化的Simpson求積symsx ifdisp('Error!Thereissomething

disp('Whendisp('Thediscretesolutionofy(x) for disp(['Theumerroris: title('Solveequationbynumericalintergration');legend('Acuratesolution','Approximatesolution');gridWhen Thediscretesolutionofy(x) Theumerroris:0.055676。08x3

x0x首先將計算定積分的Romberg方法編寫成子函數(shù)的形式functionr2=romberg梯形積分然后由梯形計算Simpson積分由Simpson積分計算Cotes積分后由ote積分計算ombeg積分，若計算的積分值滿足精度要求則推出計算，主函數(shù)中將矩陣,b,,ermbeg(,bf,)Romberg%Romberg計算定積分functionr2=romberg(a,b,f,e) whileer>e|i<5 % Simpson %Cotes Romberg symsxf disp('Thevalueoftheintergrationis:'); Thevalueoftheintergrationis:

10

)d0