《數(shù)值分析》課程設(shè)計(jì)—作業(yè)

《數(shù)值分析》課程設(shè)計(jì)—作業(yè) （本文檔內(nèi)的有些運(yùn)行結(jié)果，限于篇幅，使文檔結(jié)構(gòu)更和諧、緊湊，已做相關(guān)的改動(dòng)，程序代碼沒(méi)變） 實(shí)驗(yàn)一 1.1 水手、猴子和椰子問(wèn)題：五個(gè)水手帶了一只猴子來(lái)到南太平洋的一個(gè)荒島上，發(fā)現(xiàn)那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊，很快就入睡了。第一個(gè)水手醒來(lái)后，把椰子平分成五堆，將多余的一只給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個(gè)水手也陸續(xù)起來(lái)，和第一個(gè)水手一樣，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再給猴子，試問(wèn)原先共有幾只椰子？ 試分析椰子數(shù)目的變化規(guī)律，利用逆向遞推的方法求解這一問(wèn)題。 解： 1、 問(wèn)題分析：對(duì)于本題，比較簡(jiǎn)單，我們只需要判斷原來(lái)椰子的個(gè)數(shù)及每個(gè)人私藏了一份之后剩下的是否能被5除余1，直到最后分完。 function fentao(n) a=cat(1,7); for j=n:-1:1 a(1)=j;i=1; while i<7 a(i+1)=4*(a(i)-1)/5; i=i+1; end if a(7)==fix(a(7)) a, end end end 2、 問(wèn)題求解：通過(guò)matlab建立M文件，有如下程序： n=input(input n:); for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end if p==fix(p) break; end end disp([x,p]) 或者 input n:2000 1023 15621 >> fentao(20001) a = 15621 12496 9996 7996 6396 5116 4092 對(duì)于第一個(gè)程序，n取2000；對(duì)于第二個(gè)程序，n取20001，就能得到我們想要的結(jié)果，即原先一共有15621個(gè)椰子，最終平均每人得4092個(gè)椰子。 1.2 當(dāng)時(shí)，選擇穩(wěn)定的算法計(jì)算積分. 解： 一、問(wèn)題分析：由知： 以及： 得遞推關(guān)系： ， 但是通過(guò)仔細(xì)觀察就能知道上述遞推公式每一步都將誤差放大十倍，即使初始誤 差很小，但是誤差的傳播會(huì)逐步擴(kuò)大，也就是說(shuō)用 它構(gòu)造的算法是不穩(wěn)定的，因此我們改進(jìn)上述遞推公式（算法）如下： 通過(guò)比較不難得出該誤差是逐步縮小的，即算法是穩(wěn)定的。 2、 問(wèn)題求解：為了利用上面穩(wěn)定的算法，需要我們估計(jì)初值的值。 因?yàn)? 所以當(dāng)n=100的時(shí)，我們有： 9.090909090909091e-0049.645173882220725e-004 于是可取他們的平均值，有=9.368041486564908e-004，利用上述穩(wěn)定算法，可求得相應(yīng)的值和程序如下（限于篇幅，這里只給出了前后各十個(gè)連續(xù)的數(shù)據(jù)，詳細(xì)的數(shù)據(jù)會(huì)連同作業(yè)一并上交）： 積分計(jì)算對(duì)照表 1 n 計(jì)算值 準(zhǔn)確值 誤差 1 0.059182761 0.059182761 0 2 0.040292947 0.040292947 0 3 0.029402893 0.029402893 1.17997E-16 4 0.022708713 0.022708713 1.5278E-16 5 0.018333956 0.018333956 0 6 0.01531285 0.01531285 1.13285E-16 7 0.013125037 0.013125037 0 8 0.0114765 0.0114765 0 9 0.010193063 0.010193063 0 91 0.00099999 0.00099999 2.68937E-11 92 0.00098911 0.00098911 2.71896E-10 93 0.000978464 0.000978464 2.74854E-09 94 0.000968045 0.000968045 2.77813E-08 95 0.000957845 0.000957846 2.80771E-07 96 0.000947862 0.000947859 2.83729E-06 97 0.000938051 0.000938078 2.86687E-05 98 0.000928766 0.000928497 0.000289646 99 0.000916421 0.00091911 0.002926039 100 0.000936804 0.000909911 0.029556214 I=cat(1,100); J=cat(1,100); K=cat(1,100); I(100)=9.368041486564908e-004; format long; %求近似值 for n=99:-1:1 I(n)=((1-exp(-n))/n-I(n+1))/10; end %求精確值 for n=1:100 syms x; k=n; J(k)=int(exp(-k*x)/(exp(-x)+10),x,0,1); end %求誤差 for n=1:100 K(n)=abs((I(n)-J(n))/J(n)); end n=1:100; A=[n;I;J;K];B=A xlswrite(1_2.xls,B) >> format long >> min=(1-exp(-100))/(11*100), max=(1-exp(-100))/(100*(exp(-1)+10)) min = 9.090909090909091e-004 max = 9.645173882220725e-004 >> I100=(min+max)/2 I100 = 9.368041486564908e-004 1.3 繪制靜態(tài)和動(dòng)態(tài)的Koch分形曲線 問(wèn)題描述：從一條直線段開(kāi)始，將線段中間的三分之一部分用一個(gè)等邊三角形的另兩條邊代替，形成具有5個(gè)結(jié)點(diǎn)的新的圖形；在新的圖形中，又將圖中每一直線段中間的三分之一部分都用一個(gè)等邊三角形的另兩條邊代替，再次形成新的圖形，這時(shí)，圖形中共有17個(gè)結(jié)點(diǎn)。這種迭代繼續(xù)進(jìn)行下去可以形成Koch分形曲線。在迭代過(guò)程中，圖形中的結(jié)點(diǎn)將越來(lái)越多，而曲線最終顯示細(xì)節(jié)的多少取決于所進(jìn)行的迭代次數(shù)和顯示系統(tǒng)的分辨率。Koch分形曲線的繪制與算法設(shè)計(jì)和計(jì)算機(jī)實(shí)現(xiàn)相關(guān)。 圖1.1 Koch曲線的形成過(guò)程 解： 一、算法分析： 考慮由直線段（2個(gè)點(diǎn)）產(chǎn)生第一個(gè)圖形（5個(gè)點(diǎn)）的過(guò)程。上圖中，設(shè)和分別為原始直線段的兩個(gè)端點(diǎn)，現(xiàn)需要在直線段的中間依次插入三個(gè)點(diǎn)，，。顯然位于線段三分之一處，位于線段三分之二處，點(diǎn)的位置可看成是由點(diǎn)以點(diǎn)為軸心，逆時(shí)針旋轉(zhuǎn)600而得。旋轉(zhuǎn)由正交矩陣 實(shí)現(xiàn)。 算法根據(jù)初始數(shù)據(jù)（和點(diǎn)的坐標(biāo)），產(chǎn)生圖1中5個(gè)結(jié)點(diǎn)的坐標(biāo)。結(jié)點(diǎn)的坐標(biāo)數(shù)組形成一個(gè)矩陣，矩陣的第一行為的坐標(biāo)，第二行為的坐標(biāo)……，第五行為的坐標(biāo)。矩陣的第一列元素分別為5個(gè)結(jié)點(diǎn)的坐標(biāo)，第二列元素分別為5個(gè)結(jié)點(diǎn)的坐標(biāo)。 進(jìn)一步考慮該曲線形成過(guò)程中結(jié)點(diǎn)數(shù)目的變化規(guī)律。設(shè)第次迭代產(chǎn)生的結(jié)點(diǎn)數(shù)為，第次迭代產(chǎn)生的結(jié)點(diǎn)數(shù)為，則和中間的遞推關(guān)系為。 2、 問(wèn)題求解: ⑴、程序代碼 p=[0 0;10 0]; %P為初始兩個(gè)點(diǎn)的坐標(biāo),第一列為x坐標(biāo),第二列為y坐標(biāo) n=2; %n為結(jié)點(diǎn)數(shù) A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; %旋轉(zhuǎn)矩陣 for k=1:6 d=diff(p)/3; %diff計(jì)算相鄰兩個(gè)點(diǎn)的坐標(biāo)之差,得到相鄰兩點(diǎn)確定的向量 %則d就計(jì)算出每個(gè)向量長(zhǎng)度的三分之一,與題中將線段三等分對(duì)應(yīng) m=4*n-3; %迭代公式 q=p(1:n-1,:); %以原點(diǎn)為起點(diǎn),前n-1個(gè)點(diǎn)的坐標(biāo)為終點(diǎn)形成向量 p(5:4:m,:)=p(2:n,:); %迭代后處于4k+1位置上的點(diǎn)的坐標(biāo)為迭代前的相應(yīng)坐標(biāo) p(2:4:m,:)=q+d; %用向量方法計(jì)算迭代后處于4k+2位置上的點(diǎn)的坐標(biāo) p(3:4:m,:)=q+d+d*A; %用向量方法計(jì)算迭代后處于4k+3位置上的點(diǎn)的坐標(biāo) p(4:4:m,:)=q+2*d; %用向量方法計(jì)算迭代后處于4k位置上的點(diǎn)的坐標(biāo) n=m; subplot(2,3,k) plot(p(:,1),p(:,2),r) %繪出每相鄰兩個(gè)點(diǎn)的連線 axis([0,3,0,3]) %迭代后新的結(jié)點(diǎn)數(shù)目 %axis off end ⑵、運(yùn)行結(jié)果： 實(shí)驗(yàn)二 2.1 小行星軌道問(wèn)題：一天文學(xué)家要確定一顆小行星繞太陽(yáng)運(yùn)行的軌道，他在軌道平面內(nèi)建立以太陽(yáng)為原點(diǎn)的直角坐標(biāo)系，在五個(gè)不同的對(duì)小行星作了五次觀察，測(cè)得軌道上五個(gè)點(diǎn)的坐標(biāo)數(shù)據(jù)（單位：萬(wàn)公里）如下表所示： P1 P2 P3 P4 P5 X坐標(biāo) 53605 58460 62859 66662 68894 Y坐標(biāo) 6026 11179 16954 23492 68894 由開(kāi)普勒第一定律知，小行星軌道為一橢圓，橢圓的一般方程可表示為： 現(xiàn)需要建立橢圓的方程以供研究。 (1) 分別將五個(gè)點(diǎn)的數(shù)據(jù)代入橢圓一般方程中，寫(xiě)出五個(gè)待定系數(shù)滿(mǎn)足的等式，整理后寫(xiě)出線性方程組 以及方程組的系數(shù)矩陣和右端項(xiàng)b； (2) 用MARLAB求低階方程的指令A(yù)＼b求出待定系數(shù); (3) 分別用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系數(shù). 解： （1）、①、程序代碼： X=[53605 58460 62859 66662 68894]; Y=[6026 11197 16954 23492 68894]; B=[X.*X;2*X.*Y;Y.*Y;2*X;2*Y]; %B是系數(shù)矩陣的轉(zhuǎn)置 A=B;%A是系數(shù)矩陣 A=vpa(A,16) syms a1 a2 a3 a4 a5; a=[a1;a2;a3;a4;a5]; b=[-1;-1;-1;-1;-1]; b=sym(A*a) ②、運(yùn)行結(jié)果： A = [ 2873496025.0, 646047460.0, 36312676.0, 107210.0, 12052.0] [ 3417571600.0, 1309153240.0, 125372809.0, 116920.0, 22394.0] [ 3951253881.0, 2131422972.0, 287438116.0, 125718.0, 33908.0] [ 4443822244.0, 3132047408.0, 551874064.0, 133324.0, 46984.0] [ 4746383236.0, 9492766472.0, 4746383236.0, 137788.0, 137788.0] b = 2873496025.0*a1 + 646047460.0*a2 + 36312676.0*a3 + 107210.0*a4 + 12052.0*a5 3417571600.0*a1 + 1309153240.0*a2 + 125372809.0*a3 + 116920.0*a4 + 22394.0*a5 3951253881.0*a1 + 2131422972.0*a2 + 287438116.0*a3 + 125718.0*a4 + 33908.0*a5 4443822244.0*a1 + 3132047408.0*a2 + 551874064.0*a3 + 133324.0*a4 + 46984.0*a5 4746383236.0*a1 + 9492766472.0*a2 + 4746383236.0*a3 + 137788.0*a4 + 137788.0*a5 即是方程組的形式如下: -1=2873496025.0*a1 + 646047460.0*a2 + 36312676.0*a3 + 107210.0*a4 + 12052.0*a5 -1=3417571600.0*a1 + 1309153240.0*a2 + 125372809.0*a3 + 116920.0*a4 + 22394.0*a5 -1=3951253881.0*a1 + 2131422972.0*a2 + 287438116.0*a3 + 125718.0*a4 + 33908.0*a5 -1=4443822244.0*a1 + 3132047408.0*a2 + 551874064.0*a3 + 133324.0*a4 + 46984.0*a5 -1=4746383236.0*a1 + 9492766472.0*a2 + 4746383236.0*a3 + 137788.0*a4 + 137788.0*a5 （2） 、用指令A(yù)＼b求出待定系數(shù)的值如下： >> b=[-1;-1;-1;-1;-1];A;a=A\b a = 0.000000000080151337137572070705034479362305 -0.000000000099919136582879577308894350837601 -0.00000000017639794583181340912897589334276 -0.000012556065177690492529375516352892 0.000015497775048603393791640753306649 （3） 、1）、直接法 ①、程序代碼： function [RA,RB,n,X]=bilizy(A,b) B=[A,b];n=length(b); RA=rank(A);RB=rank(B); rank_cha=RB-RA; if rank_cha>0, disp(因?yàn)镽A!=RB，所以此方程組無(wú)解.) return end if RA==RB if RA==n disp(因?yàn)镽A=RB=n,所以此方程組有唯一解.) X=zeros(n,1);C=zeros(1,n+1); T=0;D=B;s=max(abs(D(1:n,:))); S=s; for p=1:n-1 [Y,j]=max((abs(B(p:n,p)))./S(p:n)); J=min(j); C=B(p,:);B(p,:)=B(J+p-1,:);B(J+p-1,:)=C; T=S(p);S(p)=S(J+p-1);S(J+p-1)=T; for k=p+1:n m=B(k,p)/B(p,p); B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1); A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp(請(qǐng)注意：因?yàn)镽A=RB<n,所以此方程組有無(wú)窮多解.) end end ②、運(yùn)行結(jié)果： >>b=[-1;-1;-1;-1;-1];[RA,RB,n,X]=bilizy(A,b) 因?yàn)镽A=RB=n,所以此方程組有唯一解. RA = 5 RB = 5 n = 5 X = 1.0e-004 * 0.0000 -0.0000 -0.0000 -0.1256 0.1550 >> X=vpa(X,16) X = 0.00000000008015133713756955 -0.0000000000999191365828799 -0.0000000001763979458318154 -0.00001255606517769044 0.00001549777504860351 2）、Jacobi迭代法 ①、程序代碼： %P:范數(shù)的名稱(chēng)，P=1，2，inf %error:近似解a的誤差 %maxl:迭代的最大次數(shù) function [a]=Jacobi(A,b,a0,P,error,maxl) [n n]=size(A); a=zeros(n,1); for k=1:maxl for j=1:n a(j)=(b(j)-A(j,[1:j-1,j+1:n]) a0([1:j-1,j+1:n]))/A(j,j); end erra=norm(a-a0,P); a0=a; if(erra<error) disp(近似解a分別是：) a=vpa(a,16); return end end if(erra>=error) disp(請(qǐng)注意：Jacobi迭代次數(shù)已經(jīng)超過(guò)最大迭代次數(shù)maxl。) end ②、運(yùn)行結(jié)果： >>X=[53605 58460 62859 66662 68894]; Y=[6026 11197 16954 23492 68894]; B=[X.*X;2*X.*Y;Y.*Y;2*X;2*Y]; %B是系數(shù)矩陣的轉(zhuǎn)置 A=B;%A是系數(shù)矩陣 b=[-1;-1;-1;-1;-1]; a0=[0;0;0;0;0]; >> [a]=Jacobi(A,b,a0,1,0.0001,100) 近似解a分別是： a = -0.0000000003480081375786834 -0.0000000007638525188999265 -0.000000003479009721870011 -0.000007500525036752573 -0.000007257526054518536 3）、Gauss-Seidel迭代法 ①、程序代碼： %P:范數(shù)的名稱(chēng)，P=1,2,inf;error：近似解x的誤差;maxl：迭代的最大次數(shù) function [a]=GS(A,b,a0,P,error,maxl) [n n]=size(A); a=zeros(n,1); for k=1:maxl for j=1:n aa=0; for i=1:n if i<j aa=aa+A(j,i)*a(i); end if i>j aa=aa+A(j,i)*a0(i); end end a(j)=(b(i)-aa)/A(j,j); end erra=norm(a-a0,P); a0=a; if (erra<error) disp(近似解a分別是：) a=vpa(a,16); return end end if(erra>=error) disp(請(qǐng)注意：Gauss-Seidel迭代次數(shù)已經(jīng)超過(guò)最大迭代次數(shù)maxl.) end ②、運(yùn)行結(jié)果： >> X=[53605 58460 62859 66662 68894];Y=[6026 11197 16954 23492 68894]; B=[X.*X;2*X.*Y;Y.*Y;2*X;2*Y]; A=B;b=[-1;-1;-1;-1;-1];a0=[0;0;0;0;0]; >> [a]=GS(A,b,a0,inf,0.00001,100) 近似解a分別是： a = 0.0000000001678142444990526 0.000000002024408102712696 0.000000000369101766105352 -0.000009750490974607461 -0.0001554717527254895 2.3 設(shè) (1) 將矩陣A進(jìn)行LU分解A=LU，其中U是上三角矩陣，L是主對(duì)角線上的元素都是1的下三角矩陣。 解： >> A=[20 2 3;1 8 1;2 -3 15]; >> [L,U]=lu(A) L = 1.0000 0 0 0.0500 1.0000 0 0.1000 -0.4051 1.0000 U = 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443 (2) 利用上述分解分別求解方程組 并由此求出逆矩陣。 >> b1=[1;0;0];b2=[0;1;0];b3=[0;0;1]; >> X1=U\(L\b1), X2=U\(L\b2),X3=U\(L\b3),inv(A) X1 = 0.0517 -0.0055 -0.0080 解： inv(A) = 0.0517 -0.0164 -0.0093 -0.0055 0.1237 -0.0072 -0.0080 0.0269 0.0665 X3 = -0.0093 -0.0072 0.0665 X2 = -0.0164 0.1237 0.0269 (3) 用LU分解求下列線性方程組的解 （方程組的精確解是） 解： >> A=[4 2 -3 -1 2 1 0 0 0 0; 8 6 -5 -3 6 5 0 1 0 0; 4 2 -2 -1 3 2 -1 0 3 1; 0 -2 1 5 -1 3 -1 1 9 4; -4 2 6 -1 6 7 -3 3 2 3; 8 6 -8 5 7 17 2 6 -3 5; 0 2 -1 3 -4 2 5 3 0 1; 16 10 -11 -9 17 34 2 -1 2 2; 4 6 2 -7 13 9 2 0 12 4; 0 0 -1 8 -3 -24 -8 6 3 -1]; >> b=[5;12;3;2;3;46;13;38;19;-21]; >> [L U]=lu(A); >> x=U\(L\b) x = 1.0000 -1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000 -1.0000 2.0000 結(jié)果是： 2.4 (1) 用追趕法求解方程組 (a) ， (b) ， 解： (a)、程序代碼： %追趕法 function x=zhuiganfa(A,b) n=rank(A); for i=1:n if A(i,i)==0 return; end end d=ones(n,1);a=ones(n-1,1);c=ones(n-1); for i=1:n-1 a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i); end d(n,1)=A(n,n); %解Ly=b的解，解保存在b中 for i=2:n d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1); end %求解Uy=x的解 x(n,1)=b(n,1)/d(n,1); for i=n-1:-1:1 x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1); end 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1666 0.1668 0.1662 0.1683 0.1607 0.1891 0.0829 0.4793 x = 0.4793 0.0829 0.1891 0.1607 0.1683 0.1662 0.1668 0.1666 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 運(yùn)行結(jié)果： >> clear all >> A=zeros(30,30); A(30,30)=4; b=ones(30,1); b(1,1)=2;b(30,1)=2; for i=1:29 A(i,i+1)=1;A(i+1,i)=1; A(i,i)=4; end >> x=zhuiganfa(A,b) >> A=zeros(20,20); A(19,19)=5;A(20,19)=-2; A(19,20)=-2;A(20,20)=5; b=zeros(20,1); b(1,1)=1; for i=1:18 A(i,i+2)=1;A(i+2,i)=1; A(i,i+1)=-2;A(i+1,i)=-2; A(i,i)=5; end >> [L,U]=lu(A); >> x=vpa(U\(L\b),6) 0.0000504376 0.000106309 0.0000292324 -0.0000143431 -0.0000126677 -0.00000146436 0.00000249126 0.00000131198 -9.33542*10^(-8) -2.99738*10^(-7) x = 0.242121 0.0943914 -0.0218242 -0.0313623 -0.00694256 0.00488693 0.00358612 0.000214823 -0.000784527 -0.000357862 (b)、 (2) 設(shè)計(jì)實(shí)驗(yàn)驗(yàn)證Hilbert矩陣的病態(tài)性，其中 format rat; for n=1:10 A=hilb(n); tiaojianshuA=norm(A,1)*norm(inv(A),1) end 解： tiaojianshuA = 1 tiaojianshuA = 27 tiaojianshuA = 748 tiaojianshuA = 28375 tiaojianshuA = 943656 tiaojianshuA = 29070279 tiaojianshuA = 985194890 tiaojianshuA = 33872791987 tiaojianshuA = 1099651961846 tiaojianshuA = 35353690673203 從上面的結(jié)果可以看出，Hilbert矩陣的條件數(shù)的增長(zhǎng)速度驚人，才十階，矩陣的條件數(shù)就為35353690673203，這是一個(gè)，可怕的數(shù)字，由它組成的方程組肯定是病態(tài)的。 實(shí)驗(yàn)三 3.1 用Taylor級(jí)數(shù)的第項(xiàng)多項(xiàng)式，分別在區(qū)間和上逼近正弦函數(shù)，并用計(jì)算機(jī)繪出上面31個(gè)函數(shù)的圖形。 解： 一、程序代碼： syms x; y=sin(x); ezplot(y,[0,pi]); grid on; axis([0,pi,0,1.5]); hold on; for m=1:2:30 p=taylor(y,x,m); ezplot(p,[0,pi]); axis([0,pi,0,1.5]); end 二、函數(shù)圖形： 3.2 已知四個(gè)點(diǎn)A，B，C，D的具體位置A(0, 0)，B(0, 3), C(8, 1), D(10, 5),求兩個(gè)點(diǎn)H1，H2的具體位置，使AH1+BH1+H1H2+H2C+H2D為最短。 解：function [ d ] =ju_li(a,b) d=sqrt(sum((a-b).^2)); return end syms x1 y1 x2 y2; A=[0 0];B=[0 3];C=[8 1];D=[10 5];H1=[x1 y1];H2=[x2 y2]; sum_d=ju_li(A,H1)+ju_li(B,H1)+ju_li(H1,H2)+ju_li(C,H2)+ju_li(D,H2) [x1 y1 x2 y2]= solve(diff(sum_d,x1),diff(sum_d,y1),diff(sum_d,x2),diff(sum_d,y2),x1,y1,x2,y2) sum_d = 3*x1^2 - 2*x1*x2 + 3*x2^2 - 36*x2 + 3*y1^2 - 2*y1*y2 - 6*y1 + 3*y2^2 - 12*y2 + 199 x1 = 9/4 y1 = 27/4 x2 = 15/8 y2 = 21/8 實(shí)驗(yàn)四 4.1 方程的根實(shí)際上是兩個(gè)函數(shù)的交點(diǎn)的橫坐標(biāo)，用計(jì)算機(jī)繪出兩個(gè)函數(shù)在區(qū)間的圖形，觀察圖形，分析它們的交點(diǎn)分布規(guī)律及特點(diǎn)，試寫(xiě)出方程的全部實(shí)根所在的隔根區(qū)間；并用二分法求出每一個(gè)根的近似值。 解： 一、畫(huà)函數(shù)圖像如下： >> syms x y1 y2; >> y1=sin(x);y2=1/x; >> ezplot(y1,[-6,6]);grid on; hold on; ezplot(y2,[-6,6]) >> gtext(y1=sin(x)) >> gtext(y2=1/x) 2、 用二分法求近似解: 通過(guò)上面的函數(shù)圖像可以知道，y1=sin(x)與y2=1/x在區(qū)間[-6,6]上共有四個(gè)交點(diǎn)，即是函數(shù)y=xsin(x)在區(qū)間[-6,6]上共有四個(gè)零點(diǎn)，調(diào)用matlab中的函數(shù)【>> ginput(4) 】 可以大概得到四個(gè)根的大概區(qū)間為：[-4,-2]、[-2,-0.5]、[0.5,2]、[2,4]，通過(guò)二分法求四個(gè)點(diǎn)的近似值的過(guò)程如下： %用二分法求非線性方程f(x)=0的根；fun為函數(shù)的表達(dá)式 %a,b為左右端點(diǎn)，eps為精度，x為近似根，k為二分次數(shù) function [x,k]=eff(fun,a,b,eps) if nargin<4 eps=1e-5;%如果輸入自變量數(shù)目<4，默認(rèn)eps=1e-5 end fa=feval(fun,a);fb=feval(fun,b); if fa*fb>0 disp([a,b]不是有根區(qū)間，請(qǐng)重新調(diào)整); return; end k=0; while abs(b-a)/2>eps x=(a+b)/2;fx=feval(fun,x); if fx*fa<0 b=x;fb=fx; else a=x;fa=fx; end k=k+1; end x=(a+b)/2; ①、程序代碼： ②、運(yùn)行結(jié)果： >> fun=@(x)(x*sin(x)-1) ;eps=1e-5; >> [x1,k1]=eff(fun,-4,-2,eps),[x2,k2]=eff(fun,-2,-0.5,eps), [x3,k3]=eff(fun,0.5,2,eps),[x4,k4]=eff(fun,2,4,eps) x1 = -2.7726 k1 = 17 x4 = 2.7726 k4 = 17 x3 = 1.1142 k3 = 17 x2 = -1.1142 k2 = 17 4.2 設(shè)方程有3個(gè)實(shí)根 盡可能多地采用下面六種不同計(jì)算格式，求的根或或的近似值，并觀察相應(yīng)格式的收斂性。 (1) (2) (3) (4) (5) (6) 解： >> syms x fun; >> fun=x^3-3*x-1; >> ezplot(fun,[-2.5,2.5]),grid 1、 函數(shù)圖像： 得到該函數(shù)的零點(diǎn)大概在：-1.5、-0.3、1.8的附近。 2、 迭代格式： 程序代碼: %迭代法 %用迭代法求非線性方程f(x)=0的根，fun為函數(shù)ψ（x）的表達(dá)式 %x0為初值，eps為精度（默認(rèn)1e-5），N為最大迭代次數(shù)（默認(rèn)500），x為近似根 function [x,k]=ddf(fun,x0,eps,Nmax) if nargin<4 Nmax=500; end if nargin<3 eps=1e-5; end x=x0;x0=x+2*eps;k=0; while abs(x0-x)>eps&k<Nmax x0=x;x=feval(fun,x0); k=k+1 end if k==Nmax warning(已迭代次數(shù)上限); end (1) >> fun=inline((3*x+1)/x^2); >> [x1,k1]=ddf(fun,-1.5,1e-5),[x2,k2]=ddf(fun,-0.3,1e-5),[x3,k3]=ddf(fun,1.8,1e-5), -1.5321 k1 = 28 x2 = NaN k2 = 31 x3 = NaN k3 = 48 x2 = NaN k2 = 31 x3 = NaN k3 = 48 x1 = -1.5321 k1 = 28 (2) >> fun=inline((x^3-1)/3); >> [x1,k1]=ddf(fun,-1.5,1e-5),[x2,k2]=ddf(fun,-0.3,1e-5),[x3,k3]=ddf(fun,1.8,1e-5), x2 = -0.3473 k2 = 5 x3 = -0.3473 k3 = 9 x1 = -0.3473 k1 = 12 (3) >> fun=inline((3*x+1)^(1/3)); >> [x1,k1]=ddf(fun,-1.5,1e-5),[x2,k2]=ddf(fun,-0.3,1e-5),[x3,k3]=ddf(fun,1.8,1e-5), x3 = 1.8794 k3 = 8 x2 = 1.8794 k2 = 12 x1 = 1.8794 + 0.0000i k1 = 12 >> fun=inline(1/(x^2-3)); >> [x1,k1]=ddf(fun,-1.5,1e-5),[x2,k2]=ddf(fun,-0.3,1e-5),[x3,k3]=ddf(fun,1.8,1e-5), (4) x2 = -0.3473 k2 = 5 x3 = -0.3473 k3 = 7 x1 = -0.3473 k1 = 8 (5) >> fun=inline((3+1/x)^(1/3)); >> [x1,k1]=ddf(fun,-1.5,1e-5),[x2,k2]=ddf(fun,-0.3,1e-5),[x3,k3]=ddf(fun,1.8,1e-5), x1 = 1.5396 k1 = 6 x2 = 1.5396 + 0.0000i k2 = 7 x3 = 1.5396 k3 = 5 (6) >> fun=inline(x-(1/3)*(x^3-3*x-1)/(x^2-1)); >> [x1,k1]=ddf(fun,-1.5,1e-5),[x2,k2]=ddf(fun,-0.3,1e-5),[x3,k3]=ddf(fun,1.8,1e-5), x1 = -1.5321 k1 = 3 x2 = -0.3473 k2 = 3 x3 = 1.8794 k3 = 4 4.3 一個(gè)10次多項(xiàng)式 的系數(shù)為 [1, -55, 1320, –18150, 157773, –902055, 3416930, –8409500, 12753576, –10628640, 6328800] 解： (1) 用多項(xiàng)式的求根指令roots求出方程的10個(gè)根； >> clear all >> p=[1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800]; >> roots(p) ans = 10.6051 + 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127i (3) 將方程左邊多項(xiàng)式中的9次項(xiàng)的系數(shù)–55改為–56，得到一個(gè)新的10次方程；求解新的方程，觀察根的變化是否很顯著，說(shuō)明了什么？ >> p=[1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800]; >> roots(p) ans = 21.7335 7.3501 + 7.7973i 7.3501 - 7.7973i 4.3589 + 3.3285i 4.3589 - 3.3285i 5.1831 2.4378 + 2.7974i 2.4378 - 2.7974i 0.3949 + 1.0127i 0.3949 - 1.0127i 實(shí)驗(yàn)五 5.1 給出的函數(shù)表如下： 解： ⑴、分別用Lagrange插值和Newton插值求的近似值； ①、Lagrange插值： 程序代碼： %Lagrange插值法 function f=Lagrange(x,y,x0) syms t; if length(x)==length(y) n=length(x); else disp(x和y的維數(shù)不相等！); return; end %驗(yàn)錯(cuò) f=0.0; for i=1:n l=y(i); for j=1:i-1 l=l*(t-x(j))/(x(i)-x(j)); end for j=i+1:n l=l*(t-x(j))/(x(i)-x(j)); %計(jì)算Lagrange基函數(shù) end f=f+l; %計(jì)算Lagrange插值函數(shù) simplify(f); %化簡(jiǎn) end if nargin==3 f=vpa(subs(f,t,x0),10); %計(jì)算插值點(diǎn)的函數(shù)值 else f=collect