1、PAGE 28實驗3函數(shù)式直接確定型模型PAGE 18數(shù)學實驗與MatlabPAGE PAGE 89 數(shù)學實驗與Matlab程序周曉陽華中科技大學數(shù)學系我將程序按書中的順序排列，這樣便于讀者利用。 本書程序均通過了調(diào)式?？芍苯涌截惖矫畲翱谶\行或編制M文件運行。如出現(xiàn)問題，可能是中英文標點的緣故（排版有可能使用了中文標點），請將中文標點換為英文標點試試。 實驗矩陣運算與Matlab命令1.1 知識要點與背景：知識要點和背景：日常矩陣及其運算【 A=4 2 3;1 3 2;1 3 3;3 2 2, % 表1-1、表1-2的數(shù)據(jù)分別寫成矩陣形式 B=35 20 60 45;10 15 50 40;

2、20 12 45 20 】【 C=A*B %矩陣乘法，求各訂單所對應(yīng)的原材料和勞動力 。 】 【 whos % 查看Matlab工作空間中變量及其規(guī)模 】 1.2實驗與觀察：矩陣和Matlab語言1.2.1 向量的生成和運算 【 x=linspace(0,4*pi,100); %將0,4區(qū)間100等分，產(chǎn)生了一個100維向量 y=sin(x); %計算函數(shù)值，產(chǎn)生了一個與x同維的100維函數(shù)向量y y1=sin(x).2; %計算函數(shù)向量，注意元素群運算 y2=exp(-x).*sin(x); %以x為橫坐標，y為縱坐標畫函數(shù)的圖用不同的線型將函數(shù)曲線繪制在一個圖上 plot(x,y,-,x,

3、y1,-,x,y2,.-) 】 1. 向量的創(chuàng)建 直接輸入向量。 【x1=1 2 4,x2=1,2,1,x3=x1 】 冒號創(chuàng)建向量 ?！?x1=3.4:6.7 x2=3.4:2:6.7 x3=2.6:-0.8:0 】 生成線性等分向量。 【 x=linspace(0,1,5) 】2. 向量的運算 【 y=sin(x) 】 【 y1=sin(x).2; y2=exp(-x).*sin(x); 】 1.2.2.矩陣創(chuàng)建和運算1.創(chuàng)建矩陣（1）數(shù)值矩陣的創(chuàng)建 直接輸入法創(chuàng)建簡單矩陣。【 A=1 2 3 4; 5 6 7 8; 9 10 11 12 】 【 B=-1.3,sqrt(3);(1+2)*

4、4/5,sin(5);exp(2),6 】 （2）符號矩陣的創(chuàng)建 【 syms a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34A1=a11 a12 a13 a14 ;a21 a22 a23 a24; a31 a32 a33 a34,B1=b11 b12 b13 b14 ;b21 b22 b23 b24; b31 b32 b33 b34 】 2.矩陣的運算 【 C=A1+B1,D=A1-B1 】 【 syms c cA=c*A1 】 【 C=A1*B1

5、】 ? Error using = sym/mtimes, Inner matrix dimensions must agree. 【 A2=A1(:,1:3), B1 】【 G=A2*B1 】【 g11=A2(1,:)*B1(:,1) 】 【 A, A_trans=A 】 【 H=1 2 3 ; 2 1 0 ; 1 2 3 , K=1 2 3 ; 2 1 0 ; 2 3 1h_det=det(H), k_det=det(K),H_inv=inv(H),K_inv=K-1 】 【 A=3 0 1; 1 1 0;0 1 4; B=inv(A-2*eye(3)*A, B=(A-2*eye(3)A

6、】3.分塊矩陣：矩陣的裁剪、分割、修改與抽?。?）【 A=1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1, vr=1,3;vc=1,3;A1=A(vr,vc) %取出A的1、3行和1、3列的交叉處元素構(gòu)成新矩陣A1 】 將上面的矩陣A分為四塊，并把它們賦值到矩陣B中，觀察運行后的結(jié)果?！?A11=A(1:2,1:2),A12=A(1:2,3:5),A21=A(3:4,1:2),A22=A(3:4,3:5)B=A11 A12;A21 A22 】A = 2 0 5 4 2 1 0 -1 2B = 1 2 4 -1 5 3 1 0 -1 0 2 3C = -3 4

7、18 13 13 14 20 -1 -7 -3 3 6 （2）矩陣的修改和提取 【 A=1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1 A(1,:)=0 0 0 0 0; A 】 觀察：【 B(:,2,4)= %刪除矩陣B的第2、4列 】（3）矩陣元素的抽取4.生成特殊矩陣 。 【 y1=rand(1,5), y2=rand(1,5),rand(seed,3), x1=rand(1,5), rand(seed,3), x2=rand(1,5) 】5. 常用矩陣函數(shù)6. 數(shù)據(jù)的簡單分析 【 rand(seed,1);A=rand(3,6), Asort=sort

8、(A), Amax=max(A), Asum=sum(A) 】1.2.3 Matlab工作環(huán)境和編程2.Matlab的基本設(shè)計1.3應(yīng)用、思考與練習 關(guān)系矩陣1.3.2 投入產(chǎn)出 1.3.3 循環(huán)比賽的名次 【 A=0 1 1 0; 0 0 1 1; 0 0 0 1; 1 0 0 0, e=ones(4,1); c=A*e; s=c 】 畫矩陣結(jié)構(gòu)圖的gplot指令。 （3） 【 clf, A=0 1 1 0;0 0 1 1;0 0 0 1;1 0 0 0; xy=0 1;0 0;-1 0.5;1 0.5; graphy_plot(A,xy,1,0.5), % gplot(A,xy) 】1.3

9、.4 參考程序graphy_plot.m 【 function y=graphy_plot(A,xy,l,p) %畫矩陣的有向結(jié)構(gòu)圖。 A為鄰接矩陣，xy為頂點坐標，l控制參數(shù)，l=0，畫無向圖； %l=0,畫有向圖。p為控制箭頭大小的參數(shù)。a=-max(abs(xy(:,1)*1.1;b=max(abs(xy(:,1)*1.1;c=-max(abs(xy(:,2)*1.1;d=max(abs(xy(:,2)*1.1;if l=0 gplot(A,xy),axis(a b c d),hold on,elseif l=0 U=;V=;X=;Y=; n=length(A(:,1) ; for i=

10、1:n k=find(A(i,:)=0); m=length(k); if(m=0) for j=1:m u(1)=(xy(k(j),1)-xy(i,1); v(1)=(xy(k(j),2)-xy(i,2); u(2)=eps; v(2)=eps; U=u;U;V=v;V; X=xy(i,1) xy(k(j),1);X; Y=xy(i,2) xy(k(j),2);Y; end text(xy(i,1),xy(i,2),bulletleftarrowfontsize16itV, um2str(i); hold on, end end gplot(A,xy),axis(a b c d),hold

11、on, h=quiver(X,Y,U,V,p);set(h,color,red);hold on, plot(xy(:,1),xy(:,2),k.,markersize,12),hold on,end , hold off 】實驗函數(shù)的可視化與Matlab作2.1 實驗與觀察：函數(shù)的可視化2.1.1 Matlab二維繪圖命令1.周期函數(shù)與線性p周期函數(shù) 觀察 ： 【 clf, x=linspace(0,8*pi,100);F=inline(sin(x+cos(x+sin(x);y1=sin(x+cos(x+sin(x); y2=0.2*x+sin(x+cos(x+sin(x);plot(x,y

12、1,k:,x,y2,k-) legend(sin(x+cos(x+sin(x),0.2x+sin(x+cos(x+sin(x),2) 】 2. plot指令：繪制直角坐標的二維曲線3. 圖形的屬性設(shè)置和屏幕控制 【 h=plot(0:0.1:2*pi,sin(0:0.1:2*pi); grid on set(h,LineWidth,5,color,red); set(gca,GridLineStyle,-,fontsize,16) 】 設(shè)置y坐標的刻度并加以說明，并改變字體的大小?！?h=plot(0:0.1:2*pi,sin(0:0.1:2*pi);grid on， set(gca,ytic

13、k,-1 -0.5 0 0.5 1), set(gca,yticklabel,a|b|c|d|e), set(gca,fontsize,20) 】4. 文字標注指令 【 plot(x,y1,b,x,y2,k-) , set(gca,fontsize,15,fontname,times New Roman), %設(shè)置軸對象的字體為times % New Roman，字體的大小為15 title( itPeroid and linear peroid function); %加標題，注意文字用單引號 加上 %斜杠后可輸入不同的設(shè)置，例如it表示花括號里的文字為斜體；如果有多項設(shè)置， %則可用連續(xù)輸

14、入。 xlabel(x from 0 to 8*pi itt); ylabel(ity); %說明坐標軸 text(x(49),y1(50)-0.4,fontsize15bulletleftarrowThe period function itf(x); %在坐標(x(49),y1(50)-0.4）處作文字說明， 各項設(shè)置用隔開。 %fontsize15bulletleftarrow的意義依次是：字體大小=15 畫圓點 左箭頭 text(x(14),y2(50)+1,fontsize15The linear period function itg(x)rightarrowbullet) %與上

15、一語句類似，用右箭頭 】圖2.5 文字標注 觀察指令legend和num2str的用法：在同一張圖上畫出， 這里, 并進行適當?shù)臉俗ⅰ?zxy2_2.m 【 clf, t=0:0.1:3*pi;alpha=0:0.1:3*pi; plot(t,sin(t),r-);hold on; plot(alpha,3*exp(-0.5*alpha),k:); set(gca,fontsize,15,fontname,times New Roman), xlabel(itt(deg);ylabel(itmagnitude);title( itsine wave and itAe-alphaittwave)

16、; %注意alpha的意義text(6,sin(6),fontsize15The Value itsin(t) at itt=6rightarrowbullet, HorizontalAlignment,right), %上面的語句是一整行，如果要寫成兩行，必須使用續(xù)行號 ，例如要在“ bullet,” %后換行，需寫“bullet, ”后才能換行。 % HorizontalAlignment,right 表示箭頭所指的曲線對象在 文字的右邊。text(2,3*exp(-0.5*2),fontsize15bulletleftarrow The Value of it3e-0.5 itt=,nu

17、m2str(3*exp(-0.5*2), at itt =2 ); %num2str的用法：string1 ,num2str,string2，注意方括號的使用。%legend(itsin(t),itAe-alphat) % 請結(jié)合圖形觀察此命令的使用 】 運行結(jié)果如圖2.6所示。 5. 圖形窗口的創(chuàng)建和分割 觀察：【 clf，b=2*pi;x=linspace(0,b,50);for k =1:9 y=sin(k*x); subplot(3,3,k),plot(x,y),axis(0,2*pi,-1,1)end 】x=linspace(0,4*pi,150);y1=cos(x);y=cos(c

18、os(cos(cos(cos(cos(cos(x)+sin(x)+sin(x)+sin(x);plot(x,y,x,y1,r-),grid 2.1.2多元函數(shù)的可視化與空間解析幾何（三維圖形） 本節(jié)通過高等數(shù)學的幾個例子觀察Matlab的三維繪圖功能和技巧。1. 繪制二元函數(shù) 觀察：繪制 的圖象，作定義域的裁剪。 （1）觀察meshgrid指令的效果。【 a=-0.98;b=0.98;c=-1;d=1;n=10;x=linspace(a,b,n); y=linspace(c,d,n);X,Y=meshgrid(x,y);plot(X,Y,+) 】 三維繪圖指令mesh、meshc、surf。

19、（2）做函數(shù)的定義域裁剪，觀察上述三維繪圖指令的效果。程序zxy2_4.m【 clear,clf,a=-1;b=1;c=-15;d=15;n=20;eps1=0.01;x=linspace(a,b,n);y=linspace(c,d,n);X,Y=meshgrid(x,y);for i=1:n %計算函數(shù)值z ，并作定義域裁剪 for j=1:n if (1-X(i,j)eps1|X(i,j)-Y(i,j)eps1 %if語句這樣用 z(i,j)=NaN; %作定義域裁剪，定義域以外的函數(shù)值為NaN else z(i,j)=1000*sqrt(1-X(i,j)-1.*log(X(i,j)-Y(

20、i,j); end endend zz=-20*ones(1,n);plot3(x,x,zz),grid off,hold on %畫定義域的邊界線mesh(X,Y,z) %繪圖，讀者可用meshz, surf, meshc在此替換之 xlabel(x),ylabel(y),zlabel(z), box on %把三維圖形封閉在箱體里 】 運行了zxy2_4.m 以后，有關(guān)向量存儲在工作空間中，在此基礎(chǔ)上，觀察上述等值線繪制指令的運行效果。 【 cs,h=contour(X,Y,z,15); clabel(cs,h,labelspacing,244) 】2. 三元函數(shù)可視化: slice指令

21、觀察： 繪制三元函數(shù)的可視化圖形?！?clf,x=linspace(-2,2,40); y=x; z=x; X,Y,Z=meshgrid(x,y,z); w=X.2+Y.2+Z.2;slice(X,Y,Z,w,-1,0,1,-1,0,1,-1,0,1),colorbar 】3. 空間曲線及其運動方向的表現(xiàn)：plot3和quiver指令【 clf, t=0:0.1:1.5;Vx=2*t;Vy=2*t.2;Vz=6*t.3-t.2;x=t.2;y=(2/3)*t.3;z=(6/4)*t.4-(1/3)*t.3; %由速度得到曲線plot3(x,y,z,r.-),hold on %畫飛行軌跡 %算數(shù)

22、值梯度，也就是重新計算數(shù)值速度矢量，這只是為了編程的方便，不是必須的圖2.12 飛機的飛行軌跡與方向Vx=gradient(x);Vy=gradient(y);Vz=gradient(z);quiver3(x,y,z,Vx,Vy,Vz),grid on %畫速度矢量圖xlabel(x),ylabel(y),zlabel(z) 】2.2應(yīng)用、思考和練習2.2.1 線性p周期函數(shù)zxy2_3_f.m【 function f=zxy2_3_f(x) f=sin(x+cos(x); 】zxy2_3.m【 clear，clfa=-8;b=12;n=300;xx=linspace(a,b,n);h=zxy

23、2_3_f(xx);S(1)=0;for i=2:n S(i)=S(i-1)+quad(zxy2_3_f,xx(i-1),xx(i);endsubplot(1,2,1),plot(xx,S,k-),axis(a,b,-1.5,9)subplot(1,2,2),plot(xx,h;zeros(1,length(xx),k-),axis(a,b,-1.5,1.5) 】2.2.2 平面截割法和曲面交線的繪制用平行截面法討論由曲面 構(gòu)成的馬鞍面形狀。 zxy2_6.m （ 平行截割法示例 ， 本程序的繪制兩曲面交線方法可以套用） 【 clf, a=-20;eps0=1;x,y=meshgrid(-10

24、:0.2:10); %生成平面網(wǎng)格v=-10 10 -10 10 -100 100; %設(shè)定空間坐標系的范圍colormap(gray) %將當前的顏色設(shè)置為灰色z1=(x.2-2*y.2)+eps; %計算馬鞍面函數(shù)z1=z1(x,y)z2=a*ones(size(x); %計算平面 z2=z2(x,y)r0=abs(z1-z2)=eps0; %計算一個和z1同維的函數(shù)r0,當abs(z1-z2)eps0時，r0 =0。 %可用mesh(x,y,r0)語句觀察它的圖形，體會它的作用，該方法可以套用。zz=r0.*z2;xx=r0.*x;yy=r0.*y; %計算截割的雙曲線及其對應(yīng)的坐標su

25、bplot(2,2,2), %在第2圖形窗口繪制雙曲線 h1=plot3(xx(r0=0),yy(r0=0),zz(r0=0),+); set(h1,markersize,2),hold on,axis(v),grid onsubplot(2,2,1), %在第一圖形窗口繪制馬鞍面和平面 mesh(x,y,z1);grid,hold on;mesh(x,y,z2); h2=plot3(xx(r0=0),yy(r0=0),zz(r0=0),.); %畫出二者的交線 set(h2,markersize,6),hold on,axis(v),for i=1:5 %以下程序和上面是類似的，通過循環(huán)繪制

26、一系列的平面去截割馬鞍面 a=70-i*30; %在這里改變截割平面 z2=a*ones(size(x); r0=abs(z1-z2)=1; zz=r0.*z2;yy=r0.*y;xx=r0.*x; subplot(2,2,3), mesh(x,y,z1);grid,hold on;mesh(x,y,z2);hidden off h2=plot3(xx(r0=0),yy(r0=0),zz(r0=0),.); axis(v),grid subplot(2,2,4), h4=plot3(xx(r0=0),yy(r0=0),zz(r0=0),o); set(h4,markersize,2),hold

27、 on,axis(v),grid onend 】 2.2.3 微分方程的斜率場 繪制微分方程 的斜率場，并將解曲線畫在圖中，觀察斜率場和解曲線的關(guān)系。 zxy2_5.m ( 繪制一階微分方程的斜率場和解曲線)【 clf,clear %清除當前所有圖形窗口的圖像，清除當前工作空間的內(nèi)存變量。a=0;b=4;c=0;d=4;n=15;X,Y=meshgrid(linspace(a,b,n),linspace(c,d,n); %生成區(qū)域中的網(wǎng)格。z=X.*Y; %計算斜率函數(shù)。 Fx=cos(atan(X.*Y);Fy=sqrt(1-Fx.2); %計算切線斜率矢量。quiver(X,Y,Fx,Fy

28、,0.5),hold on,axis(a,b,c,d) %在每一網(wǎng)格點畫出相應(yīng)的斜率矢量，0.5是控制矢量大小的控制參數(shù)，可以調(diào)整。x,y=ode45(zxy2_5f,0,4,0.4); %求解微分方程。 %zxy2_5f.m是方程相應(yīng)函數(shù)f(x,y)的程序，單獨編制；x0,xs=0,4為求解區(qū)間； %y0=0.4為初始值；輸出變量x,y分別為解軌線自變量和因變量數(shù)組。plot(x,y,r.-) %畫解軌線 】 zxy2_5f.m （微分方程的函數(shù)子程序）【 function dy=zxy2_5f(x,y)dy=x.*y; 】 2.2.4顏色控制和渲染及特殊繪圖指令1.地球表面的氣溫分布（sp

29、here指令） 【 a,b,c=sphere(40);t=max(max(abs(c)-abs(c);surf(a,b,c,t);axis(equal),colormap(hot), shading flat,colorbar 】2.旋轉(zhuǎn)曲面的生成：柱面指令cylinder和光照控制指令surfl 【 x=0:0.1:10;z=x;y=1./(x.3-2.*x+4); u,v,w=cylinder(y);surfl(u,v,w,45,45); shading interp 】3.若干特殊圖形 運行下面程序，了解各指令的用法和效果?！?x=1:10; y=5 6 3 4 8 1 10 3 5 6

30、;subplot(2,3,1),bar(x,y),axis(1 10 1 11)subplot(2,3,2),hist(y,x),axis(1 10 1 4)subplot(2,3,3),stem(x,y,k),axis(1 10 1 11)subplot(2,3,4),stairs(x,y,k), axis(1 10 1 11)subplot(2,3,5), x = 1 3 0.5 5;explode = 0 0 0 1;pie(x,explode)subplot(2,3,6),z=0:0.1:100; x=sin(z);y=cos(z).*10; comet3(x,y,z) 】實驗函數(shù)式直

31、接確定型模型3.1知識要點和背景：函數(shù) 直接確定性模型 3.2實驗與觀察：插值與擬合3.2.1 插值方法與多項式擬合的概念3.2.2 用Matlab作插值和擬合3.2.3 用鼠標選節(jié)點 觀察插值、擬合的效果 3.2.4 觀察程序說明zxy3_1.m【 clf,a=-1;b=1;n=100; %用內(nèi)聯(lián)函數(shù)inline命令定義函數(shù)， %在后面可直接用于函數(shù)g的計算,要改變函數(shù)做實驗，可按此格式重新定義gg=inline(x2-x4); xx=linspace(a,b,n);for i=1:n gx(i)=gxx(i); % 前面已經(jīng)用inline命令定義了g，可以這樣用g計算函數(shù)值endymin=

32、min(gx)*0.8;ymax=max(gx)*1.2;%分四個界面畫圖g的圖形，以便于結(jié)果比較subplot(2,2,1),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(近鄰插值)subplot(2,2,2),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(線性插值)subplot(2,2,3),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(樣條插值)subplot(2,2,4),plot(xx,gx,-),gri

33、d,hold on,axis(a b ymin ymax),title(多項式擬合) %用鼠標在屏幕上選點x,y,button = ginput(n)，可套用下面程序的格式button=1;x1=a;y1=gx(1);while button=1 xi,yi,button=ginput(1); subplot(2,2,1),h=plot(xi,yi,ro) %在4個圖形窗口畫點 subplot(2,2,2),h=plot(xi,yi,ro) subplot(2,2,3),h=plot(xi,yi,ro) subplot(2,2,4),h=plot(xi,yi,ro) x1=xi,x1;y1=y

34、i,y1; %將選的點存于向量x1,y1end x1=b,x1;y1=gx(n),y1; xx=linspace(a,b,n); %定義自變量xx %計算不同的插值函數(shù)：x1,y1為節(jié)點，xx為輸入自變量 ynearest=interp1(x1,y1,xx,nearest); ylinear=interp1(x1,y1,xx,linear); yspline=interp1(x1,y1,xx,spline); %多項式擬合指令p,s = polyfit(x,y,n)，n為擬合多項式次數(shù)，x,y為被擬合數(shù)據(jù)， %p為擬合多項式的系數(shù)，s是用來做誤差 估計和預測的數(shù)據(jù)結(jié)構(gòu)。 p,c=polyfit

35、(x1,y1,4); ypolyfit=polyval(p,xx); %用polyval(p,x)計算系數(shù)為p的多項式在標量或向量x處的值 subplot(2,2,1),h=plot(xx,ynearest,r-);set(h,linewidth,2) %畫圖subplot(2,2,2),h=plot(xx,ylinear,r-);set(h,linewidth,2);subplot(2,2,3),h=plot(xx,yspline,r-);set(h,linewidth,2) subplot(2,2,4),h=plot(xx,ypolyfit,r-);set(h,linewidth,2) 】

36、3.3應(yīng)用、思考和練習3.3.1若干函數(shù)的插值和擬合練習3.3.2幾個應(yīng)用問題1. 機床加工和水深流速問題 2. 內(nèi)燃機轉(zhuǎn)角與升程的關(guān)系 3. 隨高度變化的大氣壓強4. 交通事故的調(diào)查3.3.4多元函數(shù)的插值 plot(fi,x,o) fi=91 105 110 115 120 124 128; x=0 1.0869 1.9710 2.7555 3.3986 4.9073 8.3409; plot(fi,x,o) x = 0 1.0869 1.9710 2.7555 3.3986 4.9073 8.3409 1. 矩形溫箱的溫度2. 航行區(qū)域的警示線3.3.5 Fourier級數(shù)和周期函數(shù)的經(jīng)

37、驗公式 zxy3_2.m【 clf,clear,n=10;m=3;x=1:1:12;y=3.1 3.8 6.9 12.7 16.8 20.5 24.5 25.9 22.0 16.1 10.7 5.4;z=(pi/6)*x;plot(z(1:12),y(1:12),o);hold onk=1:m; %計算數(shù)據(jù)矩陣。for i=1:n X(i,:)=1 cos(k*z(i) sin(k*z(i);endc=inv(X*X)*X*y(1:n), %求解。z1=linspace(0,2*pi,30);s=; %開始計算F級數(shù)的部分和。for i=1:30; sd=1 cos(k*z1(i) sin(k

38、*z1(i); %注意k是向量。 sd=c.*sd; s=s,sum(sd);endplot(z1,s,r-),hold on, xlabel(月份),ylabel(平均氣溫) 】3.3.6由實驗到理論:從開普勒定律和牛頓萬有引力定律x=0 3 5 7 9 11 12 13 14 15;y=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x1=0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9;y1=3.1950, 3.2299, 3.2532, 3.2611, 3.2516, 3.2282, 3.1870, 3.1266,

39、 3.0594, 2.9757;clf,subplot(1,2,1),plot(x,y,k.),axis(-0.5 16 -0.5 4),subplot(1,2,2),plot(x1,y1,k.) T=0.241 0.615 1.000 1.881 11.862 29.457; R=0.387 0.723 1.000 1.524 5.203 9.539; T=T.2, R=R.3 plot(T,R,o) 運行結(jié)果（圖形略，顯示為直線）T = 0.0581 0.3782 1.0000 3.5382 140.7070 867.7148R = 0.0580 0.3779 1.0000 3.5396

40、140.8515 867.9777 實驗8.線性代數(shù)實驗 PAGE 42 實驗. 微分、積分和微分方程4.1. 知識要點和背景：微積分學基本定理 對于定義在a,b區(qū)間上的函數(shù)y=f(x),將a,b區(qū)間n等分，步長h=(b-a)/n.設(shè)分點為a=x0 x1x2x=1 1 1 ;1 1 1;1 1 1,I1=cumsum(x,1),I2=cumsum(x,2)x = 1 1 1 1 1 1 1 1 1I1 = 1 1 1 2 2 2 3 3 3I2 = 1 2 3 1 2 3 1 2 3 觀察：用累積和命令cumsum計算積分。【 x=linspace(0,pi,50); y=sin(x);T=c

41、umsum(y)*pi/(50-1);I=T(50) 】 觀察梯形公式計算的效果?！?x=linspace(0,pi,50); y=sin(x);T=trapz(x,y) 】 觀察辛普森積分公式的計算效果?！?I1=quad(sin,0,pi), I2=quad8(sin,0,pi), 】 觀察：用解常微分方程的方法計算積分?！?y0=0;x,y=ode23(myfun2_2,0,pi,y0);s=length(y);y(s) 】【 y0=0;x,y=ode45(myfun2_2,0,pi,y0);s=length(y);y(s) 】4.2.3 觀察程序及其說明 zxy4_1.m （觀察方程的

42、解曲線族，圖4.1）【 clf,clear,a=0;b=pi;c=0;d=2.2;n=20;X,Y=meshgrid(linspace(a,b,n),linspace(c,d,n);z=X.*Y; Fx=cos(atan(sin(X);Fy=sin(atan(sin(X); quiver(X,Y,Fx,Fy,0.5),hold on,axis(a,b,c,d)x,y=ode45(zxy4_1f,0,pi,0); x1,y1=ode45(zxy4_1f,0,pi,0.2);x2,y2=ode45(zxy4_1f,0,pi,0.4);x3,y3=ode45(zxy4_1f,0,pi,0.6); p

43、lot(x,y,r.-,x1,y1,r.-,x2,y2,r.-,x3,y3,r.-), xlabel(x) ,ylabel(y) 】zxy4_1f.m 【 function dy=zxy4_1f(x,y) dy=sin(x); 】. 觀察 程序zxy4_2.m （圖4.14.3）【 clf, a=0;b=4*pi;n=31; x=linspace(a,b,n); y=zeros(1,n); y(1)=0; s(1)=0; for i=1:n-1 dy=myfun2_2(x(i); y(i+1)=y(i)+dy*(x(i+1)-x(i); hh(i)=myfun2_2(x(i); s(i+1)=

44、s(i)+hh(i)*(x(i+1)-x(i); enda1=0.9*a;b1=1.1*b; % 設(shè)置坐標范圍。ymin=min(y);ymax=max(y);ymin1=ymin*0.9;ymax1=ymax*1.1; subplot(2,1,1), %在第一幅圖中畫垂線,和原函數(shù)。 fplot(1-cos(x),0,pi,r:);axis(a1 b1 ymin1 ymax1),hold on, plot(x(1) x(1),ymin ymax); subplot(2,1,2), %在第二幅圖中畫被積函數(shù)圖象。 fplot(myfun2_2,a b,-k);hold on,axis(a1 b

45、1 ymin1 ymax1) for i=1:n-1 %在第I個子區(qū)間上計算。 subplot(2,1,1), plot(x(i+1) x(i+1),ymin ymax); %在各分點畫豎線。 plot(x(i) x(i+1),y(i),y(i); %畫水平線。 h=plot(x(i+1) x(i+1),y(i),y(i+1); %畫表示增量的垂線。 set(h,linewidth,3,color,b); %設(shè)置線寬和顏色。 h1=plot(x(i) x(i+1),y(i) y(i+1),.-); %畫折線，設(shè)置圖形屬性 。 set(h1,linewidth,2,markersize,18);

46、 subplot(2,1,2), xx=x(i) x(i) x(i+1) x(i+1); yy=0 hh(i) hh(i) 0; %計算矩形頂點坐標。 patch(xx,yy,0.7 0.7 0.7); %在第二幅圖中畫矩形塊并填充顏色。 plot(x(i+1) x(i+1),ymin ymax); plot(x(1) x(1),s(i),s(i+1),r,linewidth,5); %在y軸上畫面積增量線段。 plot(x(1),y(i+1),r.,Markersize,18); subplot(2,1,1),plot(x(i+1) x(i+1),ymin ymax); h=plot(x(1

47、) x(1),s(i),s(i+1); %畫相應(yīng)的輔助線段。 set(h,linestyle,-,linewidth,5,color,red); plot(x(1),y(i+1),r.,Markersize,18); plot(x(1) x(i+1),s(i+1) y(i+1),r-) pause %暫停，n大時應(yīng)該去掉。 end 】 myfun2_2.m (選擇其它的函數(shù)進行實驗，可修改此程序)【 function y=myfun2_2(x) sin(x); %y=sin(x)./(x); 】4.3 實驗與觀察()： Matlab符號微積分簡介4.3.1 創(chuàng)建符號變量4.3.2 求符號極限l

48、imit指令 觀察 ：求下列問題的極限： 【 syms x a I1=limit(sin(x)-sin(3*x)/sin(x),x,0)I2=limit(tan(x)-tan(a)/(x-a),x,a)I3=limit(3*x-5)/(x3*sin(1/x2),x,inf) 】4.3.3 求導指令diff 1.符號求導指令diff 【 syms x y f=sym(exp(-2*x)*cos(3*x(1/2)diff(f,x)g=sym(g(x,y) %建立抽象函數(shù)。f=sym(f(x,y,g(x,y) %建立復合抽象函數(shù)。diff(f,x)diff(f,x,2) 】2.數(shù)值求導指令diff

49、【 x=linspace(0,2*pi,50);y=sin(x);dydx=diff(y)./diff(x); plot(x(1:49),dydx),grid 】4.3.4 求符號積分int 【 syms x y z I1=int(sin(x*y+z),z) I2=int(1/(3+2*x+x2),x,0,1) I3=int(1/(3+2*x+x2),x,-inf,inf) 】4.3.5 化簡、提取和代入 【 syms x at=(a+x)6-(a-x)6t_expand=expand(t)t_factor=factor(t_expand)t_simplify=simplify(t) 】 觀察

50、： 將代入表達式中求值。 【 syms a b c x y a0 y0 f=a*b+c/x*y; a=a0;b=1;c=4;x=x0;y=5; t= subs(f) 】 【 syms a b c x y a0 y0 f=a*b+c/x*y; subs(f,a,b,c,x,y,a0,1,4,x0,5) 】4.4應(yīng)用、思考和練習4.4.1 追擊問題1.追擊問題的數(shù)值模擬 2. 追蹤雷達失效的情形 3. 追線問題的解析解 【 syms y a v s0 x=sym(x(y); t=sym(t(y); %定義函數(shù)關(guān)系。 f_left=-y*diff(x,y); f_right=s0+a*t-x; %方

51、程左、右邊表達式。 r_left=diff(f_left,y), r_right=diff(f_right,y) %求導 】【 syms y d r xs=simplify(dsolve(D2x=-r*sqrt(1+Dx2)/y,x(20)=0,Dx(20)=0,y) 】【 r=0.4; y=20:-0.01:0;x=1/2*(y.(1+r)*r*20(-r)+y.(1-r)*r*20r-40*r-y.(1+r)*20(-r)+y.(1-r)*20r)/(-1+r2);shg,comet(x,y) 】 4.4.2 應(yīng)用問題1.槍支的設(shè)計【clear,clfx=0.0127,0.0254,0.0

52、381,0.0508,0.0635,0.0762,0.0889,0.1016,0.1143,0.1270,0.1397,0.1524,0.3175,0.3302,0.3429,0.3556,0.3683,0.3810,0.3937,0.4064,0.4191,0.4318,0.4445,0.4572,0.1651,0.1778,0.1905,0.2032,0.2159,0.2286,0.2413,0.2540,0.2667,0.2794,0.2921,0.3048,0.4699,0.4826,0.4953,0.5080,0.5207,0.5334,0.5461,0.5588,0.5715,0.

53、5842,0.5969,0.6096;p=0.10135,0.20064,0.27303,0.31095,0.33094,0.33991,0.34474,0.33577,0.31508,0.29578,0.27717,0.26131,0.11859,0.11238,0.10687,0.10204,0.09215,0.09308,0.08894,0.08480,0.08067,0.07722,0.07377,0.07032,0.24545,0.23097,0.21718,0.20339,0.19167,0.17995,0.16823,0.15789,0.14824,0.13927,0.13238

54、,0.12548,0.06757,0.06481,0.06205,0.05929,0.05654,0.05378,0.05102,0.04826,0.04550,0.04274,0.04067,0.03861;xx,k=sort(x);pp=p(k);plot(0,xx,0,pp,-),grid on,xlabel(槍 管 長 度 x),ylabel(壓 強 p) 】 2.天然氣井的開采量 實驗.用導數(shù)作定性分析5.1知識要點：函數(shù)作圖 用導數(shù)定性描述函數(shù)【 clf,x=linspace(-8,8,30);f=(x-3).2./(4*(x-1); plot(x,f) 】【 fplot(x-3)

55、2/(4*(x-1),-8,8) 】【 clf,x=sym(x); f=(x-3)2/(4*(x-1); ezplot(f,-8,8) , title(x-3)2/(4*(x-1),fontsize,11) , xlabel(x,fontsize,11) 】 按函數(shù)繪圖步驟繪制完整的函數(shù)圖，直接用Matlab符號演算完成必須的計算。 【 df_dx=diff(x-3)2/(4*(x-1) 】 【 sym(x);factor(df_dx) 】 【 f=inline(x-3)2/(4*(x-1);X1=-1,f(-1),X2=3,f(3) 】 求符號二階導數(shù)d2y/dx2： 【 df2_dx2=d

56、iff(x-3)2/(4*(x-1),2) 】 二階導數(shù)的因式分解：【 sym(x); factor(df2_dx2) 】【 syms x f_left=limit(x-3)2/(4*(x-1),x,1,left) f_right=limit(x-3)2/(4*(x-1),x,1,right) 】【 syms x f_minus_inf=limit(x-3)2/(4*(x-1),x,-inf,right) f_plus_inf=limit(x-3)2/(4*(x-1),x,inf,left) 】 有無斜漸近線？【 syms x, a=limit(x-3)2/(4*(x-1)/x,x,inf)

57、】【 b=limit(x-3)2/(4*(x-1)-(1/4)*x,x,inf,left) 】5.2實驗與觀察：微分方程的定性解圖示5.2.1人口增長的預測1.Malthus模型2.Logistic模型【 N=dsolve(DN=r*(1-N/Nm)*N,N(t0)=N0) 】3.微分方程解的定性分析觀察1： （1）求N = N (t)的駐點和拐點。 【 syms t dN2_dt2=diff(r*(1-N(t)/Nm)*N(t),t) , dN2_dt2=factor(dN2_dt2) 】 4.用導數(shù)作穩(wěn)定性分析 下面是繪制圖5.8的參考程序。 【 clf, N=linspace(0,300

58、,50);dN=0.3134*(1-N/250).*N;plot(N,dN),hold on,plot(0 300,0,0),plot(0,250/2,250,0,0,0,o),xlabel(N,fontsize,11),ylabel(dN,fontsize,11),text(N(32),dN(32),leftarrowitd N / d t0,相點遞增右移,fontsize,11),text(125,dN(45),itdN/dt0,相點遞減左移 rightarrow,fontsize,11);h=text(251/2,1.5,itN_m/2);set(h,fontsize,11) 】5.觀察

59、程序及其說明zxy5_1.m (繪制函數(shù)圖象，圖5.3)【 clf, x=sym(x); f=(x-3)2/(4*(x-1); g=x/4-5/4;hold on,h=line(-8 8,0,0); set(h,color,red);h=line(0 0,-8,8); set(h,color,red);line(1 1,-8 8);plot(-1 1 3,-2,0,0,o),ezplot(g,-8 8); ezplot(f,-8,8), %符號函數(shù)繪圖text(-1-0.5,-2-0.5,(-1,-2);text(1,0-0.5,(1,0); text(3,0.5,(3,0);x=1.4;te

60、xt(x,subs(f),leftarrow(x-3)2/4(x-1);x=0.6;text(x,subs(f),leftarrow(x-3)2/4(x-1);x=2.5;text(x,subs(g),leftarrow斜漸近線y=x/4-5);text(1,-2,leftarrow垂直漸近線x=1);title(x-3)2/4(x-1) 】zxy5_2.m （人口預測，圖5.5）【 global p1;clf,t1=1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960