數(shù)值分析matlab上機(jī)作業(yè)報(bào)告.docVIP

一、給定向量x0，計(jì)算初等反射陣Hk。1程序功能： 給定向量x0，計(jì)算初等反射陣Hk。2基本原理： 若的分量不全為零，則由確定的鏡面反射陣H使得；當(dāng)時(shí)，由有算法：(1)輸入x，若x為零向量，則報(bào)錯(cuò)(2)將x規(guī)范化，如果M0，則報(bào)錯(cuò)同時(shí)轉(zhuǎn)出停機(jī)否則(3)計(jì)算，如果，則(4)(5)計(jì)算(6)(7)(8)按要求輸出,結(jié)束3變量說明：x輸入的n維向量；nn維向量x的維數(shù)；MM是向量x的無窮范數(shù)，即x中絕對(duì)值最大的一項(xiàng)的絕對(duì)值；p Householder初等變換陣的系數(shù)；u Householder初等變換陣的向量Us向量x的二范數(shù)；x輸入的n維向量；nn維向量x的維數(shù)；p Householder初等變換陣的系數(shù)；u Householder初等變換陣的向量Uk數(shù)k，H*x=y，使得y的第k+1項(xiàng)到最后項(xiàng)全為零；4程序代碼：(1)function p,u=holder2(x)%HOLDER2 給定向量x0,計(jì)算Householder初等變換陣的p,u%程序功能：函數(shù)holder2給定向量x0,計(jì)算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的無窮范數(shù)，即x中絕對(duì)值最大的一項(xiàng)的絕對(duì)值；if M=0 % 如果x=0，提示出錯(cuò),程序終止; disp(Error: M=0); return;else x=x/M; % 規(guī)范化end;s=norm(x); % 求x的二范數(shù)if x(1)n %如果k值溢出，報(bào)錯(cuò)； disp(Error: kn);endH=eye(n); % 初始化H，并使H(1:k,1:k)=I；p,u=holder2(x(k:n); % 得到計(jì)算Householde初等變換陣的系數(shù)、向量U；H(k:n,k:n)=eye(n-k+1)-pu*u; % 計(jì)算H(k:n,k:n)=I-pu*u；5使用示例:情形1:X為零向量 x=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檢驗(yàn): det(H)ans = -1.0000 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檢驗(yàn): 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.2673 -0.0793 0.9604 det(H)ans = -1.0000 H*xans = 4.0000 -3.7417 0.0000 0二、設(shè)A為n階矩陣，編寫用Householder變換法對(duì)矩陣A作正交分解的程序。1程序功能：給定n階矩陣A，通過本程序用Householder變換法對(duì)矩陣A作正交分解，得出AQR2基本原理： 任一實(shí)列滿秩的mn矩陣A，可以分解成兩個(gè)矩陣的乘積，即AQR，其中Q是具有法正交列向量的mn矩陣，R是非奇異的n階上三角陣。算法：(1)輸入n階矩陣A(2)對(duì)，求Househoulder初等反射陣的。(3)計(jì)算上三角陣R，仍然存儲(chǔ)在A(4)計(jì)算正交陣Q(5)按要求輸出,結(jié)束3變量說明：A輸入的n階矩陣，同時(shí)用于存儲(chǔ)上三角陣R；n矩（方）陣A的階數(shù)；QQ是具有法正交列向量的n階矩陣；p,u向量A（k:n,k），對(duì)應(yīng)初等反射陣的,uk，jj，ii循環(huán)變量；t1計(jì)算上三角陣R的系數(shù)tj；t2計(jì)算正交矩陣Q的系數(shù)ti；4程序代碼：function Q,A=qrhh(A)%QRHH 用Householder變換法對(duì)n階矩陣A作正交分解A=QR；%程序功能：函數(shù)qrhh用Householder變換法對(duì)矩陣A作正交分解A=QR；%輸入：n階矩陣A；%輸出：Q,A。Q是具有法正交列向量的n階矩陣，% A（即R）是非奇異的n階上三角陣，仍用輸入的矩陣A存儲(chǔ)。%引用函數(shù)：% 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），對(duì)應(yīng)初等反射陣的,u for jj=k:n%計(jì)算上三角陣R（仍存貯于A） t1=dot(u,A(k:n,jj)/p;%利用向量?jī)?nèi)積求和 A(k:n,jj)=A(k:n,jj)-t1*u; end for ii=1:n%計(jì)算正交矩陣Q t2=dot(u,Q(ii,k:n)/p; %利用向量?jī)?nèi)積求和 Q(ii,k:n)=Q(ii,k:n)-t2*u; endend5使用示例:(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檢驗(yàn)： 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 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 -1.8187檢驗(yàn)： 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語(yǔ)言編寫求解此辺值問題的算法程序，采用下列三種方法，并比較三種方法的計(jì)算速度。1、用SOR迭代法求解線性方程組Au=f,用試算法確定最佳松弛因子；2、用塊 Gauss-Sediel迭代法求解線性方程組Au=f；3、（預(yù)條件）共軛斜量法。用差分代替微分，對(duì)Poisson方程進(jìn)行離散化，得到五點(diǎn)格式的線性方程組寫成矩陣形式Au=f。其中一 用SOR迭代法求解線性方程組Au=f,用試算法確定最佳松弛因子。1. 基本原理：Gauss-Seidel迭代法計(jì)算簡(jiǎn)單，但是在實(shí)際計(jì)算中，其迭代矩陣的譜半徑常接近1，因此收斂很慢。為了克服這個(gè)缺點(diǎn)，引進(jìn)一個(gè)加速因子（又稱松弛因子）對(duì)Gauss-Seidel方法進(jìn)行修正加速。假設(shè)已經(jīng)計(jì)算出第k步迭代的解（i=1，2，n），要求下一步迭代的解（i=1，2，n）。首先，用Gauss-Seidel迭代格式計(jì)算然后引入松弛因子，用松弛因子對(duì)和作一個(gè)線性組合。，i=1,2,n將二者合并成為一個(gè)統(tǒng)一的計(jì)算公式：2. 算法(1)Gauss-Seidel迭代法引入松弛因子w： 五點(diǎn)格式即為： (2)計(jì)算步驟： 第一步：給松弛因子賦初值w=1.11.8,給場(chǎng)值u和場(chǎng)源b賦初值 第二步：用不同的w進(jìn)行迭代計(jì)算。置error=0;計(jì)算在計(jì)算機(jī)上采用動(dòng)態(tài)計(jì)算形式如果|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)用）%輸入場(chǎng)初值u0、場(chǎng)源b及松弛因子w，通過五點(diǎn)差分格式進(jìn)行迭代運(yùn)算，%如果第k+1次的迭代結(jié)果與第k次的差小于精度，則可以近似認(rèn)為第k+1次的迭代%結(jié)果是精確解，然后返回迭代次數(shù)k和迭代解function u,k=SOR(u,b,w) %輸出迭代結(jié)果u，及迭代次數(shù)km=length(u); %m為u的維數(shù)h=1/(m-1); %h為步長(zhǎng)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; %計(jì)算u的修正量 u(jj,j)= u(jj,j)+du; %修正u if errorabs(du), error=abs(du); end; %統(tǒng)計(jì)誤差 end end if errork, kk=k;ww=w(i);uu=u;end %把最少迭代次數(shù)付給kk,及其w,u賦給ww,uuendt=toc %統(tǒng)計(jì)程序運(yùn)算時(shí)間 4.計(jì)算結(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.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.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.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 1.0072 1.0000 1.0000 1.0000 1.0000 contourf (uu, DisplayName, uu); figure(gcf)圖一 超松弛二 用塊Jacobi迭代法求解線性方程組Au=f。1. 基本原理：對(duì)A做自然分解A=D-L-U=D-(L+U)其中D是有A的對(duì)角線元素組成的矩陣，L是由A的對(duì)角線以下元素組成的矩陣，U是由A得對(duì)角線以上元素組成的矩陣。于是將M=D,N=L+U,代入得到Dx=(L+U)x+b任取x的初值進(jìn)行迭代2. 算法：（1）Gauss-Sediel迭代法原理五點(diǎn)差分格式： 因?yàn)锳可以寫成塊狀，即： 如果把每一條線上的節(jié)點(diǎn)看作一個(gè)組，可以把Au=f表示成塊狀求解：（2）計(jì)算步驟： 第一步：給場(chǎng)值u和場(chǎng)源b賦初值，及定義 第二步：用公式，進(jìn)行迭代計(jì)算 第三步：把第k次的u賦給ub，即ub=u；然后把第k+1次的u和ub進(jìn)行比較，看是否達(dá)到精度，如果達(dá)到精度，則輸出迭代次數(shù)k和精確解u。3. 程序k,u=kuai_GaussSeidel(n)%用塊Gauss-Sediel迭代法求解正方形域上的Poisson方程邊值問題%輸入n,對(duì)x、y軸進(jìn)行n等分；先確定場(chǎng)u的邊界及場(chǎng)源b；%用k和u分別存放迭代次數(shù)和精確解function k,u=kuai_GaussSeidel(n) %對(duì)x、y軸進(jìn)行n等分h=1/n; %步長(zhǎng)u=zeros(n+1,n+1); %對(duì)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); %計(jì)算場(chǎng)源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);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é)果的對(duì)應(yīng)元素的最大誤差 if error format short n=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.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.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 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,對(duì)x、y軸進(jìn)行n等分；先確定場(chǎng)u的邊界及場(chǎng)源b；%用k和u分別返回迭代次數(shù)和精確解function k,u=J_CG(n)h=1/n; %h步長(zhǎng)u=zeros(n+1,n+1); %對(duì)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); %計(jì)算場(chǎng)源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);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=kk/mm; %beita是 p=z+beita*p; %對(duì)p進(jìn)行更新，確定搜索方向 error=sqrt(norm(u-ub); %error是統(tǒng)計(jì)誤差 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.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 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 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四 三種迭代法效率分析 有場(chǎng)的等值圖可以看出，三種迭代方法的結(jié)果（達(dá)到相同的精度）比較一致，但是J-PCG只需要37次（耗時(shí)0.1100秒）迭代即可達(dá)到比較好的結(jié)果；超松弛迭代法需要48次（耗時(shí)0.0310秒）迭代即可達(dá)到滿意的結(jié)果；塊Gauss-Sediel迭代法需要93次（耗時(shí)1.3280）迭代才到達(dá)滿意的結(jié)果。所以對(duì)此問題來說，超松弛迭代法和J-PCG法的效率要好于塊Gauss-Sediel迭代法。一. 任務(wù)： 用MATLAB語(yǔ)言編寫連續(xù)函數(shù)最佳平方逼近的算法程序（函數(shù)式M文件）。并用此程序進(jìn)行數(shù)值試驗(yàn)，寫出計(jì)算實(shí)習(xí)報(bào)告。二. 程序功能要求： 在書中Page355或Page345的程序leastp.m（見附一）的基礎(chǔ)上進(jìn)行修改，使其更加完善。要求算法程序可以適應(yīng)不同的具體函數(shù)，具有一定的通用性。所編程序具有以下功能：1. 用Lengendre多項(xiàng)式做基，并適合于構(gòu)造任意次數(shù)的最佳平方逼近多項(xiàng)式?？衫眠f推關(guān)系 2. 被逼近函數(shù)f(x)不用內(nèi)聯(lián)函數(shù)構(gòu)造，而改用M文件建立數(shù)學(xué)函數(shù)。這樣，此程序可通過修改建立數(shù)學(xué)函數(shù)的M文件以適用不同的被逼近函數(shù)（要學(xué)會(huì)用函數(shù)句柄）。3. 要考慮一般的情況。因此，程序中要有變量代換的功能。4. 計(jì)算組合系數(shù)時(shí)，計(jì)算函數(shù)的積分采用變步長(zhǎng)復(fù)化梯形求積法（見附三）。5. 程序中應(yīng)包括幫助文本和必要的注釋語(yǔ)句。另外，程序中也要有必要的反饋信息。6. 程序輸入：（1）待求的被逼近函數(shù)值的數(shù)據(jù)點(diǎn)(可以是一個(gè)數(shù)值或向量)（2）區(qū)間端點(diǎn)：a,b。7. 程序輸出：（1）擬合系數(shù)：（2）待求的被逼近函數(shù)值1. 基本原理：設(shè)內(nèi)積空間，取線性無關(guān)組構(gòu)造的子空間 span,對(duì)被逼近函數(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在中的最佳平方逼近元。若取正交函數(shù)組對(duì)權(quán)函數(shù)(x)滿足=dx=（i,j=0,1,，n）以Legendre正交多項(xiàng)式為函數(shù)的最佳平方逼近多項(xiàng)式，若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)的三次最佳平方逼近多項(xiàng)式 :s(x)=c1p1(x)+c2 p2(x)+c3 p3(x)+c4 p4(x)2 程序功能：該程序可對(duì)被擬合函數(shù)做以Legendre多項(xiàng)式為基的最佳平方擬合，擬合數(shù)據(jù)可以是一個(gè)點(diǎn)或一個(gè)向量，并繪出該函數(shù)與擬合多項(xiàng)式以及泰勒展式在指定區(qū)間上的對(duì)比圖。3 使用說明：程序用法簡(jiǎn)單，只需輸入擬合數(shù)據(jù)點(diǎn)x，擬合區(qū)間a、b的值，以及c,s=leasp(x,a,b),程序即可自動(dòng)運(yùn)行。4 源程序：(1) 主程序M文件%LEASTP.M: 以Legendre多項(xiàng)式為基的最佳平方擬合function c,s=leastp(x,a,b)%x:輸入擬合點(diǎn)%func:輸入擬合函數(shù)%c:輸出擬合系數(shù)%s:輸出擬合值if nargin3, disp(必須輸入三個(gè)參數(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)*(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)*p