matlab上機(jī)作業(yè)報告(計算初等反射陣-用Householder變換法對矩陣A作正交分解-連續(xù)函數(shù)最佳平方逼近等)

1、一、給定向量x0，計算初等反射陣Hk。1程序功能： 給定向量x0，計算初等反射陣Hk。2基本原理： 若的分量不全為零，則由確定的鏡面反射陣H使得；當(dāng)時，由有算法：(1)輸入x，若x為零向量，則報錯(2)將x規(guī)范化，如果M0，則報錯同時轉(zhuǎn)出停機(jī)否則(3)計算，如果，則(4)(5)計算(6)(7)(8)按要求輸出,結(jié)束3變量說明：x輸入的n維向量；nn維向量x的維數(shù)；MM是向量x的無窮范數(shù)，即x中絕對值最大的一項的絕對值；p Householder初等變換陣的系數(shù)；u Householder初等變換陣的向量Us向量x的二范數(shù)；x輸入的n維向量；nn維向量x的維數(shù)；p Householder初等變換

2、陣的系數(shù)；u Householder初等變換陣的向量Uk數(shù)k，H*x=y，使得y的第k+1項到最后項全為零；4程序代碼：(1)function p,u=holder2(x)%HOLDER2 給定向量x0,計算Householder初等變換陣的p,u%程序功能：函數(shù)holder2給定向量x0,計算Householder初等變換陣的p,u；%輸入：n維向量x；%輸出：p,u。p是Householder初等變換陣的系數(shù)，% u是Householder初等變換陣的向量U。n=length(x); % 得到n維向量x的維數(shù)；p=1;u=0; % 初始化p,u；M=max(abs(x); % 得到向量x的

3、無窮范數(shù)，即x中絕對值最大的一項的絕對值；if M=0 % 如果x=0，提示出錯,程序終止; disp(Error: M=0); return;else x=x/M; % 規(guī)范化end;s=norm(x); % 求x的二范數(shù)if x(1)n %如果k值溢出，報錯； disp(Error: kn);endH=eye(n); % 初始化H，并使H(1:k,1:k)=I；p,u=holder2(x(k:n); % 得到計算Householde初等變換陣的系數(shù)、向量U；H(k:n,k:n)=eye(n-k+1)-pu*u; % 計算H(k:n,k:n)=I-pu*u；5使用示例:情形1:X為零向量 x

4、=0,0,0,0; H=holderk(x,1)Error: M=0H = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1情形2：K值溢出： x=1,2,3,4; H=holderk(x,5)Error: kn情形3：K值為1： x=2,3,4,5; H=holderk(x,1)H = -0.2722 -0.4082 -0.5443 -0.6804 -0.4082 0.8690 -0.1747 -0.2184 -0.5443 -0.1747 0.7671 -0.2911 -0.6804 -0.2184 -0.2911 0.6361檢驗: det(H)ans = -1.0000

5、H*xans = -7.3485 0.0000 0.0000 0.0000情形4：(1)K值為3： x=4,3,2,1; H=holderk(x,3)H = 1.0000 0 0 0 0 1.0000 0 0 0 0 -0.8944 -0.4472 0 0 -0.4472 0.8944檢驗: det(H)ans = -1 H*xans = 4.0000 3.0000 -2.2361 0(2)K值為2: x=4,3,2,1; H=holderk(x,2)H = 1.0000 0 0 0 0 -0.8018 -0.5345 -0.2673 0 -0.5345 0.8414 -0.0793 0 -0

6、.2673 -0.0793 0.9604 det(H)ans = -1.0000 H*xans = 4.0000 -3.7417 0.0000 0二、設(shè)A為n階矩陣，編寫用Householder變換法對矩陣A作正交分解的程序。1程序功能：給定n階矩陣A，通過本程序用Householder變換法對矩陣A作正交分解，得出AQR2基本原理： 任一實列滿秩的mn矩陣A，可以分解成兩個矩陣的乘積，即AQR，其中Q是具有法正交列向量的mn矩陣，R是非奇異的n階上三角陣。算法：(1)輸入n階矩陣A(2)對，求Househoulder初等反射陣的。(3)計算上三角陣R，仍然存儲在A(4)計算正交陣Q(5)按要

7、求輸出,結(jié)束3變量說明：A輸入的n階矩陣，同時用于存儲上三角陣R；n矩（方）陣A的階數(shù)；QQ是具有法正交列向量的n階矩陣；p,u向量A（k:n,k），對應(yīng)初等反射陣的,uk，jj，ii循環(huán)變量；t1計算上三角陣R的系數(shù)tj；t2計算正交矩陣Q的系數(shù)ti；4程序代碼：function Q,A=qrhh(A)%QRHH 用Householder變換法對n階矩陣A作正交分解A=QR；%程序功能：函數(shù)qrhh用Householder變換法對矩陣A作正交分解A=QR；%輸入：n階矩陣A；%輸出：Q,A。Q是具有法正交列向量的n階矩陣，% A（即R）是非奇異的n階上三角陣，仍用輸入的矩陣A存儲。%引用函數(shù)

8、：% holder2；示例 p,u=holder2(x);n,n=size(A);%求矩（方）陣A的階數(shù)；Q=eye(n);%構(gòu)造正交矩陣Q(1)=I；for k=1:n-1 p,u=holder2(A(k:n,k);%向量A（k:n,k），對應(yīng)初等反射陣的,u for jj=k:n%計算上三角陣R（仍存貯于A） t1=dot(u,A(k:n,jj)/p;%利用向量內(nèi)積求和 A(k:n,jj)=A(k:n,jj)-t1*u; end for ii=1:n%計算正交矩陣Q t2=dot(u,Q(ii,k:n)/p; %利用向量內(nèi)積求和 Q(ii,k:n)=Q(ii,k:n)-t2*u; ende

9、nd5使用示例:(1)A為3階矩陣： A=1 2 3; 2 3 0; 3 4 5A = 1 2 3 2 3 0 3 4 5 q,r=qrhh(A)q = -0.2673 0.8729 0.4082 -0.5345 0.2182 -0.8165 -0.8018 -0.4364 0.4082r = -3.7417 -5.3452 -4.8107 0 0.6547 0.4364 -0.0000 0.0000 3.2660檢驗： q*rans = 1.0000 2.0000 3.0000 2.0000 3.0000 0.00003.0000 4.0000 5.0000(2)A為4階矩陣： A=1 2

10、3 4; 2 3 0 1; 3 4 5 6;1 6 8 0A = 1 2 3 4 2 3 0 1 3 4 5 6 1 6 8 0 q,r=qrhh(A)q = -0.2582 0.0597 -0.2660 -0.9268 -0.5164 -0.1045 0.8434 -0.1049 -0.7746 -0.2688 -0.4662 0.3323 -0.2582 0.9556 -0.0222 0.1399r = -3.8730 -6.7132 -6.7132 -6.1968 0 4.4647 6.4805 -1.4783 0 -0.0000 -3.3070 -3.0178 0 0.0000 0 -

11、1.8187檢驗： q*rans = 1.0000 2.0000 3.0000 4.0000 2.0000 3.0000 -0.0000 1.0000 3.0000 4.0000 5.0000 6.0000 1.0000 6.0000 8.0000 0數(shù)值求解正方形域上的Poisson方程邊值問題用MATLAB語言編寫求解此辺值問題的算法程序，采用下列三種方法，并比較三種方法的計算速度。1、用SOR迭代法求解線性方程組Au=f,用試算法確定最佳松弛因子；2、用塊 Gauss-Sediel迭代法求解線性方程組Au=f；3、（預(yù)條件）共軛斜量法。用差分代替微分，對Poisson方程進(jìn)行離散化，得到

12、五點格式的線性方程組寫成矩陣形式Au=f。其中一 用SOR迭代法求解線性方程組Au=f,用試算法確定最佳松弛因子。1. 基本原理：Gauss-Seidel迭代法計算簡單，但是在實際計算中，其迭代矩陣的譜半徑常接近1，因此收斂很慢。為了克服這個缺點，引進(jìn)一個加速因子（又稱松弛因子）對Gauss-Seidel方法進(jìn)行修正加速。假設(shè)已經(jīng)計算出第k步迭代的解（i=1，2，n），要求下一步迭代的解（i=1，2，n）。首先，用Gauss-Seidel迭代格式計算然后引入松弛因子，用松弛因子對和作一個線性組合。，i=1,2,n將二者合并成為一個統(tǒng)一的計算公式：算法(1)Gauss-Seidel迭代法引入松弛

13、因子w： 五點格式即為： (2)計算步驟： 第一步：給松弛因子賦初值w=1.11.8,給場值u和場源b賦初值 第二步：用不同的w進(jìn)行迭代計算。置error=0;計算在計算機(jī)上采用動態(tài)計算形式如果|du|error則error=|du|,如果errore,則停機(jī)，輸出結(jié)果u,k.第三步：比較不同的w的迭代次數(shù)，用kk存放最小迭代次數(shù)，用ww和uu存放相應(yīng)的w及u。3 程序 u,k=SOR(u,b,w) %（被下面程序調(diào)用）%輸入場初值u0、場源b及松弛因子w，通過五點差分格式進(jìn)行迭代運算，%如果第k+1次的迭代結(jié)果與第k次的差小于精度，則可以近似認(rèn)為第k+1次的迭代%結(jié)果是精確解，然后返回迭代次

14、數(shù)k和迭代解function u,k=SOR(u,b,w) %輸出迭代結(jié)果u，及迭代次數(shù)km=length(u); %m為u的維數(shù)h=1/(m-1); %h為步長N=10000;e=0.0000001; %e為精度for k=1:N %k為記錄迭代次數(shù) error=0; for j=2:m-1 for jj=2:m-1 sum=4*u(jj,j)-u(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1); du=w*(h2*b(jj,j)-sum)/4; %計算u的修正量 u(jj,j)= u(jj,j)+du; %修正u if errorabs(du), error=a

15、bs(du); end; %統(tǒng)計誤差 end end if errork, kk=k;ww=w(i);uu=u;end %把最少迭代次數(shù)付給kk,及其w,u賦給ww,uuendt=toc %統(tǒng)計程序運算時間 4.計算結(jié)果： format short n=10; kk,ww,uu=SOR_5dianchafen(n) t = 0.0310kk = 48ww = 1.6000uu = Columns 1 through 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.00

16、39 1.0044 1.0047 1.0045 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0

17、203 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.

18、0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000 contourf (uu, DisplayName, uu); figure(gcf)圖一 超松弛二 用塊Jacobi迭代法求解線性方程組Au=f。 基本原理：對A做自然分解A=D-L-U=D-(L+U)其中D是有A的對角線元素組成的矩陣，

19、L是由A的對角線以下元素組成的矩陣，U是由A得對角線以上元素組成的矩陣。于是將M=D,N=L+U,代入得到Dx=(L+U)x+b任取x的初值進(jìn)行迭代 算法：（1）Gauss-Sediel迭代法原理五點差分格式： 因為A可以寫成塊狀，即： 如果把每一條線上的節(jié)點看作一個組，可以把Au=f表示成塊狀求解：（2）計算步驟： 第一步：給場值u和場源b賦初值，及定義 第二步：用公式，進(jìn)行迭代計算 第三步：把第k次的u賦給ub，即ub=u；然后把第k+1次的u和ub進(jìn)行比較，看是否達(dá)到精度，如果達(dá)到精度，則輸出迭代次數(shù)k和精確解u。程序k,u=kuai_GaussSeidel(n)%用塊Gauss-Sed

20、iel迭代法求解正方形域上的Poisson方程邊值問題%輸入n,對x、y軸進(jìn)行n等分；先確定場u的邊界及場源b；%用k和u分別存放迭代次數(shù)和精確解function k,u=kuai_GaussSeidel(n) %對x、y軸進(jìn)行n等分h=1/n; %步長u=zeros(n+1,n+1); %對u賦初值u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;b=zeros(n+1,n+1); %計算場源bfor i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h2; endend b=h2*b;for i=2:

21、n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1);b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);endA=zeros(n-1,n-1); %定義矩陣的子塊Afor i=1:n-1 if i1, A(i,i-1)=-1; end if i2&jn, u(2:n,j)=pinv(A)*(u(2:n,j-1)+u(2:n,j+1)+b(2:n,j);end end error=max(max(abs(u-ub); %error是前后兩次迭代結(jié)果的對應(yīng)元素的最大誤差 if error format short n=

22、10; k,u=kuai_GaussSeidel(n)t = 1.3280k = 93u = Columns 1 through 8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045 1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0000 1.

23、0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0000 1.0023 1.0045 1.0066 1

24、.0084 1.0100 1.0110 1.0113 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103

25、1.0072 1.0000 1.0000 1.0000 1.0000 contourf (u, DisplayName, u); figure(gcf)圖二 塊Gauss-Sediel迭代法三 （預(yù)條件）共軛斜量法求解線性方程組Au=f。1.基本原理（1）預(yù)條件共軛斜量法原理（2）預(yù)優(yōu)矩陣的選取 2. 算法3.程序%用J-PCG求解正方形域上的Poisson方程邊值問題%輸入n,對x、y軸進(jìn)行n等分；先確定場u的邊界及場源b；%用k和u分別返回迭代次數(shù)和精確解function k,u=J_CG(n)h=1/n; %h步長u=zeros(n+1,n+1); %對u賦初值u(1,1:n+1)=1;

26、u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;b=zeros(n+1,n+1); %計算場源bfor i=2:n+1 for j=2:n+1 b(i,j)=(i-1)*(j-1)*h2; endend b=h2*b;for i=2:n b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1); b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);endz=zeros(n-1,n-1);r=zeros(n-1,n+1);p=zeros(n-1,n-1);bb=zeros(n-1,n-1)

27、;A=zeros(n-1,n-1); %定義矩陣的子塊Afor i=1:n-1 if i1, A(i,i-1)=-1; end if i2&j1&i1&in-1,r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i+1)-p(:,i-1);end if i=n-1,r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i-1);end end z(:,1:n-1)=pinv(A)*r(:,1:n-1); % kk=0;mm=0; for i=1:n-1 kk=kk+dot(z(:,i),r(:,i); mm=mm+dot(zz(:,i),rr(:,i); end beita

28、=kk/mm; %beita是 p=z+beita*p; %對p進(jìn)行更新，確定搜索方向 error=sqrt(norm(u-ub); %error是統(tǒng)計誤差 if error format short n=10; k,u=J_CG(n)t = 0.1100k = 37u = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045 1.0037 1.0023 1.0000 1.0000 1

29、.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088 1.0073 1.0045 1.0000 1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128 1.0107 1.0066 1.0000 1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162 1.0136 1.0084 1.0000 1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189 1.0160 1.0100 1.0000

30、1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203 1.0174 1.0110 1.0000 1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201 1.0175 1.0113 1.0000 1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175 1.0155 1.0103 1.0000 1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113 1.0103 1.0072

31、 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 contourf (u, DisplayName, u); figure(gcf) 圖三 J-PCG四 三種迭代法效率分析 有場的等值圖可以看出，三種迭代方法的結(jié)果（達(dá)到相同的精度）比較一致，但是J-PCG只需要37次（耗時0.1100秒）迭代即可達(dá)到比較好的結(jié)果；超松弛迭代法需要48次（耗時0.0310秒）迭代即可達(dá)到滿意的結(jié)果；塊Gauss-Sediel迭代法需要93次（耗時1.3280）迭代才到達(dá)滿意的結(jié)果。所以對此問題

32、來說，超松弛迭代法和J-PCG法的效率要好于塊Gauss-Sediel迭代法。一. 任務(wù)： 用MATLAB語言編寫連續(xù)函數(shù)最佳平方逼近的算法程序（函數(shù)式M文件）。并用此程序進(jìn)行數(shù)值試驗，寫出計算實習(xí)報告。二. 程序功能要求： 在書中Page355或Page345的程序leastp.m（見附一）的基礎(chǔ)上進(jìn)行修改，使其更加完善。要求算法程序可以適應(yīng)不同的具體函數(shù)，具有一定的通用性。所編程序具有以下功能：用Lengendre多項式做基，并適合于構(gòu)造任意次數(shù)的最佳平方逼近多項式?？衫眠f推關(guān)系 被逼近函數(shù)f(x)不用內(nèi)聯(lián)函數(shù)構(gòu)造，而改用M文件建立數(shù)學(xué)函數(shù)。這樣，此程序可通過修改建立數(shù)學(xué)函數(shù)的M文件以適

33、用不同的被逼近函數(shù)（要學(xué)會用函數(shù)句柄）。要考慮一般的情況。因此，程序中要有變量代換的功能。計算組合系數(shù)時，計算函數(shù)的積分采用變步長復(fù)化梯形求積法（見附三）。程序中應(yīng)包括幫助文本和必要的注釋語句。另外，程序中也要有必要的反饋信息。程序輸入：（1）待求的被逼近函數(shù)值的數(shù)據(jù)點(可以是一個數(shù)值或向量)（2）區(qū)間端點：a,b。7. 程序輸出：（1）擬合系數(shù)：（2）待求的被逼近函數(shù)值1. 基本原理：設(shè)內(nèi)積空間，取線性無關(guān)組構(gòu)造的子空間 span,對被逼近函數(shù)f(x)不屬于構(gòu)造逼近函數(shù)s(x):s(x)= 求出系數(shù)c0,c1, ,cn,使f(x)-s(x)22=min,則函數(shù)s(x)就是f(x)H在中的最佳

34、平方逼近元。若取正交函數(shù)組對權(quán)函數(shù)(x)滿足=dx=（i,j=0,1,，n）以Legendre正交多項式為函數(shù)的最佳平方逼近多項式，若f(x)=xcos(x),x-1，1，用p1(x)=1,p2(x)=x,p3(x)=(3x2-1)/2,p4(x)=(5x3-3x)/2構(gòu)造函數(shù)f(x)的三次最佳平方逼近多項式 :s(x)=c1p1(x)+c2 p2(x)+c3 p3(x)+c4 p4(x)2 程序功能：該程序可對被擬合函數(shù)做以Legendre多項式為基的最佳平方擬合，擬合數(shù)據(jù)可以是一個點或一個向量，并繪出該函數(shù)與擬合多項式以及泰勒展式在指定區(qū)間上的對比圖。3 使用說明：程序用法簡單，只需輸入擬

35、合數(shù)據(jù)點x，擬合區(qū)間a、b的值，以及c,s=leasp(x,a,b),程序即可自動運行。4 源程序：(1) 主程序M文件%LEASTP.M: 以Legendre多項式為基的最佳平方擬合function c,s=leastp(x,a,b)%x:輸入擬合點%func:輸入擬合函數(shù)%c:輸出擬合系數(shù)%s:輸出擬合值if nargin3, disp(必須輸入三個參數(shù)); return;endff=func;for i=1:4 pp(i)=2/(2*(i-1)+1);endc(1)=quad(func,a,b)/pp(1)*(2/(b-a);c(2)=quad(fp2,a,b, , ,a,b)/pp(2

36、)*(2/(b-a);c(3)=quad(fp3,a,b, , ,a,b)/pp(3)*(2/(b-a);c(4)=quad(fp4,a,b, , ,a,b)/pp(4)*(2/(b-a);s=c(1)+c(2)*p2(x,a,b)+c(3)*p3(x,a,b)+c(4)*p4(x,a,b);x0=a:0.15:b;s0=c(1)+c(2)*p2(x0,a,b)+c(3)*p3(x0,a,b)+c(4)*p4(x0,a,b);plot(x0,s0,r-); %繪制多項式擬合圖hold onfplot(ff,a,b,feval(ff,a),feval(ff,b),c); % 繪制擬合函數(shù)圖sym

37、s xxff1=xx*cos(xx);y1=taylor(ff1,xx,4);hold onezplot(y1,a,b); % 繪制擬合函數(shù)的泰勒展式圖像title(comparison figure of fitted funcyion,fitting polynomials and taylor formulation);plot(x,s,r *); %數(shù)據(jù)點legend(原 函 數(shù) f(x),最佳逼近s(x),泰勒展開T(x),數(shù)據(jù)點);(2) legendre多項式M文件function y=p2(x,a,b)t=(2*x-a-b)./(b-a);y=t;function y=p3(x

38、,a,b)t=(2*x-a-b)/(b-a);y=(3*t.2-1)/2;function y=p4(x,a,b)t=(2*x-a-b)./(b-a);y=(5*t.3-3*t)/2;(3) 函數(shù)與多項式內(nèi)積M文件function y=fp2(x,a,b)y=feval(p2,x,a,b).*feval(func,x);function y=fp3(x,a,b)y=feval(p3,x,a,b).*feval(func,x);function y=fp4(x,a,b)y=feval(p4,x,a,b).*feval(func,x);(4) 函數(shù)M文件function f=func(x)f=x.

39、*cos(x);5 計算結(jié)果（1） X為一個數(shù)據(jù)點時c,s=squfit(0.6,-1,1)c = 0.0000 0.7174 0 -0.1821s = 0.3671輸出圖形：（2） X為一個向量時c,s=leastp(-2:0.2:2,-2,2)c = 0.0000 0.1155 0 -1.0781s = Columns 1 through 6 0.9625 0.4054 -0.0062 -0.2884 -0.4574 -0.5294 Columns 7 through 12 -0.5205 -0.4470 -0.3250 -0.1706 0.0000 0.1706 Columns 13 t

40、hrough 18 0.3250 0.4470 0.5205 0.5294 0.4574 0.2884 Columns 19 through 21 0.0062 -0.4054 -0.9625輸出圖形：6 分析比較：從圖中可看出用Legendre多項式做擬合對函數(shù)在不同的區(qū)間上都有較好的逼近，而泰勒展式只在-1，1上比較精確，而在其他區(qū)間上誤差會變得較大。附錄資料：不需要的可以自行刪除Pascal/C/C+語句對比（補(bǔ)充版）一、Hello world 先看三種語言的樣例：Pascalbegin writeln(Hello world);end.C#include int main() prin

41、tf(Hello world!n); return 0;C+#include using namespace std;int main()cout Hello world! endl; return 0; 從這三個程序可以看到一些最基本的東西。在Pascal中的begin和end，在C/C+里就是；Pascal主程序沒有返回值，而C/C+返回0（好像在C中可以為NULL）。在C/C+中，main函數(shù)以前的是頭文件，樣例中C為stdio.h，C+除了iostream還有第二行的using namespace std，這個是打開命名空間的，NOIP不會考這個，可以不管，只要知道就行了。 此外說明

42、注釋單行用/，段落的話Pascal為，C/C+為/* */。* 常用頭文件（模板）#include #include #include #include #include #include using namespace std;int main() system(“pause”);return 0;二、數(shù)據(jù)類型及定義 這里只列出常用的類型。1、整型PascalC/C+范圍shortint-128 127integershort-32768 32767longintInt -2147483648 2147483647int64long long-9223372036854775808 9223

43、372036854775807byte-0 255wordunsigned short0 65535longwordunsigned int0 4294967295qwordunsigned long long0 18446744073709551615 * 當(dāng)對long long 變量賦值時，后要加LLLong long x=6327844632743269843LL* 如果位移 x2LL* Linux: printf(“%lldn”,x);* Windows: printf(“%I64dn”,x);2、實型PascalC/C+范圍realfloat2.9E-39 1.7E38single-

44、1.5E-45 3.4E38doubledouble5.0E-324 1.7E3083、字符即字符串 字符在三種語言中都為char，C里沒有字符串，只有用字符數(shù)組來代替字符串，Pascal和C+均為string。Pascal中字符串長度有限制，為255，C+則沒有。 字符串和字符在Pascal中均用單引號注明，在C/C+中字符用單引號，字符串用雙引號。4、布爾類型 Pascal 中為 boolean，C/C+ 為 bool。值均為True 或 False。C/C+中除0外bool都為真。5、定義 常量的定義均為 const，只是在C/C+中必須要注明常量的類型。在C/C+中還可以用宏來定義常量

45、，此時不注明類型。PascalC/C+const a = 60; b = -a + 30; d = ;const int a = 60;const int b = - a + 30;const string d = “”;define MAXN 501 /這個是宏 * 宏定義其實就是直接在程序相應(yīng)的位置替換： #define randomize srand(unsigned time(NULL) #define wait for(int w=0;w a;cout a;cout a endl;特別說明C+中cin一個字符的話會自動跳過空格和回車，Pascal和C則會讀入空格和回車。在Pascal

46、中writeln(a:n:m) 表示在n個字符寬的輸出域上輸出a保留m位小數(shù)。例如：pascal write(a:6) c/c+ printf(“%6d”,a) Pascal write(a:6:2) c/c+ printf(“%6.2f”,a) C+ 如果用 cout ? (繁瑣！) 需要加頭文件 #inlude cout setprecision(2)a; /作用永久 cout setw(6)a; /作用臨時 以下三個進(jìn)制設(shè)定都是永久作用： cout deca; 相當(dāng) printf(“%d”,a); /十進(jìn)制 cout hexa; 相當(dāng) printf(“%X”,a); /十六進(jìn)制 cout

47、 octa; 相當(dāng) printf(“%o”,a); /八進(jìn)制例如：cout 12hex12oct1212endl;輸出：12c1414 C 的輸入輸出里面的字符串中%表示變量，%后面的字目表示變量類型。下面是類型表：%hd1個short型整數(shù)%d1個int型整數(shù)%u1個unsigned int型整數(shù)%I64d1個long long型整數(shù)%c1個字符%s1個C字符串%f1個float型實數(shù)%lf1個double型實數(shù)%10.4f輸出1個總寬度為10,保留4位小數(shù)的實數(shù) 文件輸入輸出：Pascalassign(input, test.in);assign(output, test.out);res

48、et(input);rewrite(output);read(a, b);writeln(a, b);close(input);close(output);CFILE *fin = fopen(“test.in”, “r”);FILE *fout = fopen(“test.out”, “w”);fscanf(fin, “%d%d”, &a, &b);fprintf(fout, “%d%d”, a, b);fclose(fin); fclose(fout);C+#include using namespace std;ifstream fin(“test.in”);ofstream fout(

49、“test.out”);fin a b;fout a b endl;fin.close(); fout.close();因為C+的讀入較慢，個人建議C+的話使用C的輸入方式。當(dāng)然也有人用C的讀入，C+的輸出的，這種方式我們稱之為城鄉(xiāng)結(jié)合。*中國計算機(jī)學(xué)會競賽須知發(fā)布的C讀寫程序：（C+ 也能用，cin,cout,scanf,printf 可混用）#include int main() int a,b; freopen(“sum.in”,”r”,stdin);freopen(“sum.out”,”w”,stdout); scanf(“%d%d”,&a,&b); printf(“%dn”,a+b)

50、; return 0; 或者：freopen(“sum.in”,”r”,stdin);freopen(“sum.out”,”w”,stdout);ios:sync_with_stdio(false); 取消同步，cin,cout的速度就不慢了！ cinab;couta+bendl; return 0;以下擴(kuò)充c/c+混用是可行的：#include #include using namespace std;int main() int a,b,c,d; freopen(sum.in,r,stdin); freopen(sum.out,w,stdout); scanf(%d%d,&a,&b); c

51、incd; printf(%dn,a+b); couta+b+c+dsn).Cwhile(scanf(%s%d,s,&n)!=EOF).四、賦值語句及運算符號 一一對應(yīng)的關(guān)系PascalC/C+賦值運算賦值:=基本運算加+減-乘*除（實數(shù)）/ (double)除法取整div(int) / (int)取余mod%比較等于=不等于!=大于大于等于=小于小于等于=邏輯且and&或or|非not!位運算左移(*2)shl且and&或or|非not異或xor其他增一inc(x)x+減一dec(x)x- 在C/C+中對某個變量自身進(jìn)行運算可以簡寫為 變量名 運算符號= 改變量 如 x += 8 就表示 x

52、 = x + 8， 即 inc(x, 8)。 在 C/C+里還存在一種三目運算 變量名 = 條件 ? 值A(chǔ) : 值B 如 x = x 0 ? x : -x; /表示若x 0 則取 x， 否則取 x， 同 if x 0 then x := x else x := -x;五、條件語句1、if C/C+中if 語句的條件必須要用括號括起來，后面不使用then。PascalC/C+if a b then flag := true else flag := false;if (a b) flag = true;else flag = false;2、多種分支 C/C+中為switch，Pascal為ca

53、se：PascalC/C+case x of 1: inc(x); 2: dec(x); else x := x * x;end;switch (x) case 1: x +; break; case 2: x -; break; default: x *= x; 切記C/C+中一定要寫break，后果你可以去掉break，運行看看就知道了。六、循環(huán)語句1、forPascalC/C+for 變量名 := 初始值 to(downto) 終止值 dofor (變量名=初始值;條件;改變方式)for i := 5 to 10 do dec(a);/終止值大于初始值用 tofor i := 5 dow

54、nto 1 do dec(a);/終止值小于于初始值用 downtofor (i = 5; i = 1; i-) a-;/*只要i 滿足條件就會一直循環(huán)。C/C+中i是實數(shù)、指針都可以*/C/C+中for的特殊用法：/變量為實數(shù)for (double i = 1; i 符號為間接引用，后面會提到。for (type1 *p = head - next; p; p = p - next) printf(“%d”, p - k);2、whilePascalC/C+while 條件 dowhile (條件)while i 0 do dec(i);while (i != 0) i-;/也可寫作 whi

55、le (i) i-;/在C/C+中非0即為真。3、repeat-until & do-whilePascalC/C+repeat 語句 until 結(jié)束條件;do while (運行條件)repeat int(i) until i 100;do i+; while (i = 100);七、數(shù)組 Pascal中數(shù)組的下標(biāo)可以隨意定義，而C/C+下標(biāo)始終為從0開始到(數(shù)組大小1)。PascalC/C+定義a : array 1.100 of integer;b :array 1.10,1.10 of int64;int a100;int b1010;含義a 為大小為100的integer數(shù)組，合法

56、下標(biāo)為1到100b 為大小為10*10的int64數(shù)組，合法下標(biāo)為1,1到10,10a 為大小為100的int數(shù)組，合法下標(biāo)為0到99b 為大小為10*10的int數(shù)組，合法下標(biāo)為0,0到9,9;使用inc(a21);b2,2:=b1,1+b1,2+b2,1;a21+;b11=b01+b00+b10; 數(shù)組清零PascalC/C+Fillchar(a, sizeof(a), 0);memset(a, 0, sizeof(a);/頭文件包含 string.h*如果要填最大： memset(a,127,sizeof(a) （但達(dá)不到 INT_MAX） 如果要填最?。?memset(a,128,si

57、zeof(a) (但達(dá)不到 INT_MIN） 如果填0： memset(a,0,sizeof(a) 如果填-1： memset(a,-1,sizeof(a)八、字符串 C風(fēng)格的字符串就是字符數(shù)組。 C+和Pascal的字符串使用基本相同，只是C+中字符串下標(biāo)以0開始，Pascal以1開始。字符串處理很多這里不一一列舉，只寫最常用的幾個。PascalC (包含)定義用：char sC+（包含）定義用：string s輸入輸出Readln(s);Writeln(s);Scanf(“%s”,s);Printf(“%sn”,s);注：不能輸入輸出c+的字符串Cins;Couts = s 的區(qū)別： ge

58、tline(cin,s)cins一次性整行讀入，直至行末尾。只讀入一個“單詞”，遇空格和行末停止。例如輸入；How are you?s=” How are you?”讀入整串含空格例如輸入；How are you?s=”How”如果三個都讀：cins1s2s3*C+ 數(shù)字與數(shù)值之間的轉(zhuǎn)換：#include #include #include /必須加入using namespace std;int main() string text = 152; int number; stringstream ss;ss number; /string - int coutnumber+100endl;

59、ss string string str = ss.str(); return 0;九、過程和函數(shù)1、過程 在C/C+中沒有過程，但可以把返回值為“空”的函數(shù)理解為過程。PascalC/C+無參過程procedure 過程名;說明部分begin 語句部分 end;/說明部分、begin、end語句部分統(tǒng)稱為過程體void 函數(shù)名(); 主體部分; return ;帶參過程procedure 過程名(形參表)過程體void 函數(shù)名(形參表)過程體 值傳和址傳：當(dāng)一個參數(shù)是值傳時，形參在子過程中相當(dāng)于一個局部變量，對它的改變不影響實在的參數(shù)值。址傳則會影響。下例中a為值傳，b為址傳。初始a = 5

60、,b = 5，運行后a = 5,b = 10；PascalC/C+var a, b:integer;procedure doit(a:integer; var b:integer);begin b := a + b; a := a + b;end;begina := 5;b := 5;doit(a, b);writeln(a, , b);end.void doit(int a, int &b) HYPERLINK a a認(rèn)為值參，b認(rèn)為變量傳參 b += a; a += b; return ;int main()int a = 5, b = 5;doit(a, b);cout a b;retu