微分、積分和微分方程
4.1. 知識要點和背景：微積分學(xué)基本定理
??????
4.2 實驗與觀察（Ⅰ）：數(shù)值微積分
4.2.1實驗：積分定義、微分方程和微積分基本定理的聯(lián)系
????????? ◆步驟1.
?????? 【?? n=4;n1=n+1;x=linspace(0,pi,n1);f=myfun2_2(x);[0:n;x;f]?? 】
?????? ◆步驟2.
?【????? y(1)=0; y(2)=y(1)+f(1)*(x(2)-x(1));
?? P_intial=[x(1),y(1)],P_final=[x(2),y(2)],???? 】
??????? ◆步驟3.
【????? y(3)=y(2)+f(2)*(x(3)-x(2));
?P_intial=[x(2),y(3)],? P_final=[x(3),y(3)]???????? 】
?????? ◆ 步驟4。
【?? Dy=y(3)-y(2)?? 】
【??? DS=f(2)*(x(3)-x(2))??? 】??
??????? ◆步驟5 .
4.2.2? 求解數(shù)值積分的Matlab積分命令
??????? ◆觀察cumsum指令的功能。
【??? x=[1 1 1 1 1]；I=cumsum(x)????? 】
??????? ◆觀察：用累積和命令cumsum計算積分。
【??? x=linspace(0,pi,50); y=sin(x);
T=cumsum(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?? （觀察方程的解曲線族，圖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);?
plot(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.1～4.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)=s(i)+hh(i)*(x(i+1)-x(i));
end
a1=0.9*a;b1=1.1*b;??????? % 設(shè)置坐標(biāo)范圍。
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 b1 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);
??? subplot(2,1,2),
?????? xx=[x(i) x(i) x(i+1) x(i+1)];? yy=[0 hh(i) hh(i) 0]; %計算矩形頂點坐標(biāo)。??????? 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) 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 求符號極限limit指令
??????? ◆觀察 ：求下列問題的極限：
?????? 【?? 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)/(x^3*sin(1/x^2))',x,inf)???? 】
4.3.3 求導(dǎo)指令diff
?1.符號求導(dǎo)指令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))')??????????????? %建立復(fù)合抽象函數(shù)。
diff(f,x)
diff(f,x,2)????????? 】
2.數(shù)值求導(dǎo)指令diff
??????? 【???? 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+x^2),x,0,1)
? I3=int(1/(3+2*x+x^2),x,-inf,inf)????????? 】
4.3.5 化簡、提取和代入
????? 【?? syms x a
t=(a+x)^6-(a-x)^6
t_expand=expand(t)
t_factor=factor(t_expand)
t_simplify=simplify(t)????????????? 】
????? ◆ 觀察： 將 代入表達式 中求值。
?【???? 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)用、思考和練習(xí)
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;????????? %方程左、右邊表達式。
?r_left=diff(f_left,y), r_right=diff(f_right,y)? %求導(dǎo)??????????????????? 】
【???? syms y d r
?xs=simplify(dsolve('D2x=-r*sqrt(1+Dx^2)/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*20^r-40*r-y.^(1+r)*20^(-r)+y.^(1-r)*20^r)/(-1+r^2);
shg,comet(x,y)?? 】
?
4.4.2 應(yīng)用問題
1.槍支的設(shè)計
【clear,clf
x=[0.0127,0.0254,0.0381,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.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,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.天然氣井的開采量

200多個MATLAB經(jīng)典教程和MATLAB論文請查看：matlab教程