1、上機報告1. 共軛梯度法（1）計算過程如下：第一步：取初始向量, 計算 第步：計算（2）syms n; %定義一個變量nA=input('請輸入矩陣A') b=input('請輸入矩陣b')n=size(A)X=zeros(n);D=zeros(n);R=zeros(n);x=X(:,1); %將矩陣X的第一列賦值給x作為初始向量r=b-A*x; %將x代入求得初始非零殘向量R(:,1)=r;d=r; %求得初始方向向量D(:,1)=d;for i=1:n; %利用循環(huán)進行迭代求得各向量 AF(i)=R(:,i)'*R(:,i)/(D(:,i)'

2、*A*D(:,i); X(:,i+1)= X(:,i)+AF(i)*D(:,i); R(:,i+1)=b-A*X(:,i+1); BT(i)=norm(R(:,i+1)2/norm(R(:,i)2; D(:,i+1)=R(:,i+1)+BT(i)*D(:,i) ; end x=X(:,i+1) %輸出計算結(jié)果2. 龍貝格積分s=input('請輸入被積函數(shù)表達式：f(x)=','s') %鍵盤輸入被積函數(shù)表達式 f=inline(s)a=input('請輸入積分下限a=')%輸入積分下限a b=input('請輸入積分上限b=')

3、%輸入積分上限b f0=f(a);%求下限值f(a) f1=f(b);%求上限值f(b) j=zeros(8,4);%定義一個矩陣用來存放T,S,C,R值 j(1,1)=(b-a)*0.5*(f0+f1);%計算出T1for i=2:8;%對i循環(huán)賦值 t=2(i-2); fj=zeros(t,1);%構(gòu)造一個矩陣為接下來存放f(a+(2i-1)*(b-a)/2(k+1)的值做準備 for m=1:t; fj(m,1)=f(a+(2*m-1)*(b-a)/2(i-1);%將計算得到的f(a+(2i-1)*(b-a)/2(k+1)的值賦值給對應的矩陣元素 ff=sum(fj);%對矩陣所有元素求

4、和 end j(i,1)=0.5*j(i-1,1)+(1/2(i-1)*ff;%得到所有的T值 j(i,2)=j(i,1)+(1/3)*(j(i,1)-j(i-1,1);%得到所有的S值 j(i,3)=j(i,2)+(1/3)*(j(i,2)-j(i-1,2);%得到所有的C值 j(i,4)=j(i,3)+(1/3)*(j(i,3)-j(i-1,3);%得到所有的R值 s=j(i,4)-j(i-1,4); end j %輸出計算結(jié)果表 If=vpa(j(i,4),7) %得到精確的積分結(jié)果3. 三樣條插值X=zeros(1,6) ; %定義一些下面將用到的矩陣y=zeros(1,6);f=ze

5、ros(6);l=zeros(1,6);for i=1:1:6 X(i)=(i-1)*2)/5-1; %將區(qū)間等分取點 y(i)=1/(1+25*X(i)*X(i); %得到對應點的函數(shù)值endfor j=2:6 f(1,1)= y(1); f(j, 1)=y(j);for k=2:j f(j,k)=(f(j,k-1)-f(j-1,k-1)/(X(j)-X(j-k+1); %利用循環(huán)求差商endendf %得到差商表，其中對角線上的為需要的差商值syms x;l=x,x,x,x,x,x;g=l-X;n=g;for t=2:6 n(1,t)= n(1,t-1)* n(1,t);endN=zero

6、s(1,1);N=f(1,1);for r=1:5N= N +n(1,r)*f(r+1,r+1); %利用循環(huán)求5次牛頓插值多項式endN %得到5次牛頓插值多項式P=zeros(1,101);X1=-1:0.02:1; %取間距為0.02的等分點作為作圖的橫坐標for u=1:101 x=X1(u);subs(N) ;p(1,u)= subs(N(1,1) ; %用上面所求出的5次牛頓插值多項式算得所取等分點的函數(shù)值作為縱坐標endplot(X1,p, '-c') %繪制5次牛頓插值曲線hold on XX=zeros(1,11) ; %定義一些下面將用到的矩陣yy=zero

7、s(1,11);F=zeros(11);L=zeros(1,11);for I=1:1:11 XX(I)=(I-1)*2)/10-1; %將區(qū)間等分取點 yy(I)=1/(1+25*XX(I)*XX(I); %得到對應點的函數(shù)值endfor J=2:11 F(1,1)= yy(1); F(J, 1)=yy(J);for K=2:J F(J,K)=(F(J,K-1)-F(J-1,K-1)/(XX(J)-XX(J-K+1); %利用循環(huán)求差商endendF %得到差商表，其中對角線上的為需要用的差商值syms x1;L=x1,x1,x1,x1,x1,x1,x1,x1,x1,x1,x1;G=L-XX

8、;M=G;for T=2:11 M(1,T)= M(1,T-1)* M(1,T);endH=F(1,1);for R=1:10H= H +M(1,R)*F(R+1,R+1); %利用循環(huán)求10次牛頓插值多項式endH %得到10次牛頓插值多項式P=zeros(1,101);XX1=-1:0.02:1; %取間距為0.02的等分點作為作圖的橫坐標for u1=1:101 x1=XX1(u1);subs(H) ;P(1,u1)= subs(H) ; %用上面所求出的10次牛頓插值多項式算得所取等分點的函數(shù)值作為縱坐標endplot(XX1,P, '-g')hold on %繪制10

9、次牛頓插值曲線x3=-1:0.02:1; %取間距為0.02的等分點作為作圖的橫坐標y3=1./(1+25*x3.2); %用原函數(shù)計算得到所取等分點的縱坐標plot(x3,y3, '-k')hold on %繪制原函數(shù)f(x)的曲線h=0.2;u=0.5;z=0.5; d=zeros(11,1);d(1)=6*(yy(2)-yy(1)/0.2-50/(26*26)/0.2;d(11)=6*(-50/(26*26)-(yy(11)-yy(10)/0.2)/0.2;a=0.5*ones(1,10);b=2*ones(1,11);c=diag(a,1)+diag(a,-1)+dia

10、g(b); c(1,2)=1;c(11,10)=1;e=zeros(11,1); %得到三次樣條插值第二種邊界條件對應的嚴格三對角占優(yōu)矩陣e(1,1)=d(1);e(11,1)=d(11);for ii=2:10 e(ii,1)=F(ii+1,3); ende; %得到d0,d1,.dn用矩陣e表示ff=horzcat(c,e) %利用d0,d1,.dn與上面所得的三對角占優(yōu)矩陣構(gòu)造增廣矩陣ffy0=zeros(11,1); %定義一些下面將要用的矩陣uu=zeros(11,12);uu(1,1)=ff(1,1);ll=zeros(11,11);ll(1,1)=1;y0(1,1)=ff(1,1

11、2);for tt=2:11ll(tt,tt)=1;ll(tt,1)=ff(tt,1)/ uu(1,1);l(tt,tt-1)=ff(tt,tt-1) /uu(tt-1,tt-1); for jj=tt:12uu(1,jj)= ff(1,jj);uu(tt,jj)=ff(tt,jj)- ll(tt,tt-1)*uu(tt-1,jj); %對增廣矩陣進行LU分解，同時進行追法求y0y0(tt,1)=uu(tt,12); endendxx=zeros(11,1);for kk=10:-1:1 xx(11,1)=y0(11,1)/uu(11,12); xx(kk,1)=(y0(kk,1)-uu(kk

12、,kk+1)*xx(kk+1,1)/uu(kk,kk); %用趕法求得xxendM0=xx;Sx=zeros(10,1);for i0=1:10syms x2;Sx=sym(Sx);Sx(i0)= (XX(i0+1)-x2)*(XX(i0+1)-x2)*(XX(i0+1)-x2)/(6*h)*M0(i0)+.(x2-XX(i0)*(x2-XX(i0)*(x2-XX(i0)/(6*h)*M0(i0+1)+ .(yy(i0)-(h*h/6)*M0(i0)*(XX(i0+1)-x2)/h+.(yy(i0+1)-(h*h/6)*M0(i0+1)*(x2-XX(i0)/h; %利用循環(huán)求不同區(qū)間段的三樣

13、條插值多項式endSx %得到不同區(qū)間段的三樣條插值多項式，用矩陣表示w1=XX(1):0.01:XX(2); %將第一段等分曲線段等分成20段，得各等分點橫坐標P1=zeros(1,21); for t1=1:21 x2=w1(t1);subs(Sx(1);P1(1,t1)= subs(Sx(1); %求得第一段等分曲線段等分成20段的各等分點橫坐標對應的縱坐標endplot(w1,P1) %繪制利用三樣條插值多項式所得的第一等分段曲線hold on w2=XX(2):0.01:XX(3); %將第2段等分曲線段等分成20段，得各等分點橫坐標P2=zeros(1,21);for t2=1:2

14、1 x2=w2(t2); subs(Sx(2);P2(1,t2)= subs(Sx(2); %求得第2段等分曲線段等分成20段的各等分點橫坐標對應的縱坐標endplot(w2,P2) %繪制利用三樣條插值多項式所得的第2等分段曲線hold onw3=XX(3):0.01:XX(4); %將第3段等分曲線段等分成20段，得各等分點橫坐標P3=zeros(1,21);for t3=1:21 x2=w3(t3); subs(Sx(3);P3(1,t3)= subs(Sx(3); %求得第3段等分曲線段等分成20段的各等分點橫坐標對應的縱坐標endplot(w3,P3) %繪制利用三樣條插值多項式所得

15、的第3等分段曲線hold on w4=XX(4):0.01:XX(5); %將第4段等分曲線段等分成20段，得各等分點橫坐標P4=zeros(1,21);for t4=1:21 x2=w4(t4); subs(Sx(4);P4(1,t4)= subs(Sx(4); %求得第4段等分曲線段等分成20段的各等分點橫坐標對應的縱坐標endplot(w4,P4) %繪制利用三樣條插值多項式所得的第4等分段曲線hold onw5=XX(5):0.005:XX(6); %將第5段等分曲線段等分成40段，得各等分點橫坐標P5=zeros(1,41);for t5=1:41 x2=w5(t5); subs(S

16、x(5);P5(1,t5)= subs(Sx(5); %求得第5段等分曲線段等分成40段的各等分點橫坐標對應的縱坐標endplot(w5,P5) %繪制利用三樣條插值多項式所得的第5等分段曲線hold on%以下依次繪制其余各等分段的曲線w6=XX(6):0.005:XX(7); P6=zeros(1,41);for t6=1:41 x2=w6(t6); subs(Sx(6);P6(1,t6)= subs(Sx(6);endplot(w6,P6) %利用三樣條插值多項式繪制第6等分段曲線hold on w7=XX(7):0.01:XX(8);P7=zeros(1, 21);for t7=1:2

17、1 x2=w7(t7); subs(Sx(7);P7(1,t7)= subs(Sx(7);endplot(w7,P7) %利用三樣條插值多項式繪制第7等分段曲線hold on w8=XX(8):0.01:XX(9);P8=zeros(1,21);for t8=1:21 x2=w8(t8); subs(Sx(8);P8(1,t8)= subs(Sx(8);endplot(w8,P8) %利用三樣條插值多項式繪制第8等分段曲線hold onw9=XX(9):0.01:XX(10);P9=zeros(1,21);for t9=1:21 x2=w9(t9); subs(Sx(9);P9(1,t9)=

18、subs(Sx(9);endplot(w9,P9) %利用三樣條插值多項式繪制第9等分段曲線hold onw10=XX(10):0.01:XX(11);P10=zeros(1,21);for t10=1:21 x2=w10(t10); subs(Sx(10);P10(1,t10)= subs(Sx(10);endplot(w10,P10) %利用三樣條插值多項式繪制第10等分段曲線4. 龍格庫塔法%龍格庫塔法求3階常微分方程的初值問題的通用程序h=input('請輸入步長值h=')%通過鍵盤輸入步長ha=input('請輸入給定的x的初始點值a=')%通過鍵盤輸

19、入起始點的x值b=input('請輸入要求近似值的點的x值b=')%通過鍵盤輸入所求近似值點的x值n=(b-a)/h%求得步長個數(shù)s=input('請輸入m階微分方程轉(zhuǎn)化為一階微分方程的后的表達式f(x,y1,y2,y3)=','s')%輸入3階常微分方程轉(zhuǎn)化為一階常微分方程后的表達式f=inline(s)y=zeros(3,n+1);%定義一個矩陣來存放循環(huán)求得的y1,y2,y3值k1=zeros(3,n);%定義一個矩陣來存放循環(huán)求得的k1值k2=zeros(3,n);%定義一個矩陣來存放循環(huán)求得的k2值k3=zeros(3,n);%定義一個

20、矩陣來存放循環(huán)求得的k3值k4=zeros(3,n);%定義一個矩陣來存放循環(huán)求得的k4值A(chǔ)=input('請將微分方程在初始點a的不同階的的值按階數(shù)由低到高輸入到矩陣中A中構(gòu)成一個3x1的列矩陣A=')%將給定的初值輸入到列矩陣A中y(:,1)=A;%將y1,y2,y3的初值賦值給矩陣y的第一列for i=2:n+1%引入循環(huán) k1(1,i-1)=h*y(2,i-1);%寫出k1矩陣第一行的數(shù)學表達式k1(2,i-1)=h*y(3,i-1);%寫出k1矩陣第二行的數(shù)學表達式k1(3,i-1)=h*f(a,y(1,i-1),y(2,i-1),y(3,i-1);%寫出k1矩陣第三行的數(shù)學表達式k2(1,i-1)=h*(y(2,i-1)+k1(2,i-1)/2);%寫出k2矩陣第一行的數(shù)學表達式k2(2,i-1)=h*(y(3,i-1)+k1(3,i-1)/2);%寫出k2矩陣第二行的數(shù)學表達式k2(3,i-1)=h*f(a+(i-2)*h+h/2),y(1,i-1)+k1(1,i-1)/2,y(2,i-1)+k1(2,i-1)/2,y(3,i-1)+k1(3,i-1)/2);%寫出k2矩陣第三行的數(shù)學表達式k3(1,i-1)=h*(y(