項目四無窮級數(shù)與微分方程

程實驗1無窮級數(shù)(基礎實驗)實驗目的觀察無窮級數(shù)部分和的變化趨勢,進一步理解級數(shù)的審斂法以及冪級數(shù)部分和對逼近.掌握用Mathematica求無窮級數(shù)的和,求冪級數(shù)的收斂域,展開函數(shù)為冪級數(shù)以及展開周期函數(shù)為傅里葉級數(shù)的方法.數(shù)項級數(shù)例(教材例n=1nn=1輸入s[n_]=Sum[1/k^2,{k,n}];data=Table[s[n],{n,100}];ListPlot[data];N[Sum[1/k^2,{k,Infinity}]]N[Sum[1/k^2,{k,Infinity}],40]則輸出(1)中級數(shù)部分和的變化趨勢圖.圖級數(shù)的近似值為.輸入s[n_]=Sum[1/k,{k,n}];data=Table[s[n],{n,50}];ListPlot[data,PlotStyle->PointSize[]];則輸出(2)中級數(shù)部分和的的變化趨勢圖.圖nn=1輸入命令Clear[sn,g];sn=0;n=1;g={};m=3;While[1/n>10^-m,sn=sn+(-1)^(n-1)/n;g=Append[g,Graphics[{RGBColor[Abs[Sin[n]],0,1/n],Line[{{sn,0},{sn,1}}]}]];n++];Show[g,PlotRange->{,},Axes->True];則輸出所給級數(shù)部分和的圖形(圖),從圖中可觀察到它收斂于附近的一個數(shù).圖例求1的值.輸入Sum[x^(3k),{k,1,Infinity}]得到和函數(shù)例(教材輸入例設a=10n,nn!nClear[a];a[n_]=10^n/(n!);vals=Table[a[n],{n,1,25}];ListPlot[vals,PlotStyle->PointSize[]]則輸出a的散點圖(),從圖中可觀察a的變化趨勢.輸入Sum[a[nSum[a[n],{n,l,Infinity}]則輸出所求級數(shù)的和.圖求冪級數(shù)的收斂域n=0輸入Clear[a];a[n_]=4^(2n)*(x-3)^n/(n+1);時發(fā)散.為了求出收斂區(qū)間的端點,輸入ydd=Solve[steptwo==1,x]zdd=Solve[steptwo==-1,x]則輸出由此可知,當47<x<49時,級數(shù)收斂,當x<47或x>49時,級數(shù)發(fā)散.1616為了判斷端點的斂散性,輸入Simplify[a[n]/.x->(49/16)]則輸出右端點處冪級數(shù)的一般項為因此,在端點x=49處,級數(shù)發(fā)散.再輸入Simplify[a[n]/.x->(47/16)]則輸出左端點處冪級數(shù)的一般項為因此,在端點x=47處,級數(shù)收斂.也可以在收斂域內(nèi)求得這個級數(shù)的和函數(shù).輸入Sum[4^(2n)*(x-3)^n/(n+1),{n,0,Infinity}]則輸出函數(shù)的冪級數(shù)展開cosx輸入Series[Cos[x],{x,0,6}]則輸出注:這是帶皮亞諾余項的麥克勞林展開式.lnx在x=1處的6階泰勒展開式.輸入Series[Log[x],{x,1,6}]則輸出輸入serl=Series[ArcTan[x],{x,0,5}];Poly=Normal[serl]通過作圖把arctanx和它的近似多項式進行比較.輸入Plot[Evaluate[{ArcTan[x],Poly}],{x,-3/2,3/2},PlotStyle->{Dashing[{}],GrayLevel[0]},AspectRatio->l]則輸出所作圖形(圖),圖中虛線為函數(shù)arctanx,實線為它的近似多項式.圖多項式.輸入Clear[f];poly2=Normal[Series[f[x],{x,1,8}]]Plot[Evaluate[{f[x],poly2}],{x,,},PlotRange->{-2,3/2},PlotStyle->{Dashing[{}],GrayLevel[0]}]則得到近似多項式和它們的圖.圖xx似多項式,并形成動畫進一步觀察.所以輸入Do[Plot[{Sum[(-1)^j*x^(2j+1)/(2j+1)!,{j,0,k}],Sin[x]},{x,-40,40},PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}],{k,1,45}]則輸出為sinx的3階和91階泰勒展開的圖形.選中其中一幅圖形,雙擊后形成動畫.圖是最后一幅圖.圖例利用冪級數(shù)展開式計算5240(精確到1010).根據(jù)(1+x)m在x=0處的展開式有故前n(n>2)項部分和為輸入命令s[n_]=3(1-1/(5*3^4)-Sum[Product[5i-1,{i,1,k-1}]/(5^kk!3^(4k)),{k,2,n-1}]);r[n_]=Product[5i-1,{i,1,n-1}]/5^n/n!3^(4n-5)/80;delta=10^(-10);n0=100;Do[Print["n=",n,",","s[n]=",N[s[n],20]];If[r[n]<delta,Break[]];If[n==n0,Print["failed"]],{n,n0}]則輸出結(jié)果為傅里葉級數(shù)將g(x)展開成傅里葉級數(shù).輸入Clear[g];g[x_]:=-1/;-Pi<=x<0g[x_]:=1/;0<=x<Pig[x_]:=g[x-2Pi]/;Pi<=xPlot[g[x],{x,-Pi,5Pi},PlotStyle->{RGBColor[0,1,0]}];則輸出g(x)的圖形(圖.圖因為g(x)是奇函數(shù),所以它的傅里葉展開式中只含正弦項.輸入b2[n_]:=b2[n]=2Integrate[1*Sin[n*x],{x,0,Pi}]/Pi;fourier2[n_,x_]:=Sum[b2[k]*Sin[k*x],{k,1,n}];tu[n_]:=Plot[{g[x],Evaluate[fourier2[n,x]]},{x,-Pi,5Pi},PlotStyle->{RGBColor[0,1,0],RGBColor[1,,]},DisplayFunction->Identity];(*tu[n]是以n為參數(shù)的作圖命令*)tu2=Table[tu[n],{n,1,30,5}];toshow=Partition[tu2,2];Show[GraphicsArray[toshow]](*GraphicsArray是把圖形排列的命令*)則輸出6個排列著的圖形(圖),每兩個圖形排成一行.可以看到n越大,g(x)的傅里葉級數(shù)的前n項和與g(x)越接近.圖實驗2微分方程(基礎實驗)實驗目的理解常微分方程解的概念以及積分曲線和方向場的概念,掌握利用Mathematica求微分方程及方程組解的常用命令和方法.例(教材例求微分方程y,+2xy=xex2的通解.Clear[x,y];DSolve[y'[x]+2x*y[x]==x*Exp[-x^2],y[x],x]或DSolve[D[y[x],x]+2x*y[x]==x*Exp[-x^2],y[x],x]其中C[1]是任意常數(shù).例(教材例求微分方程xy+yex=0在初始條件y=2e下的特解.x=1輸入Clear[x,y];DSolve[{x*y'[x]+y[x]-Exp[x]==0,y[1]==2E},y[x],x]DSolve[y''[x]-2y'[x]+5y[x]==Exp[x]*Cos[2x],y[x],x]y=2x+ex輸入g1=Table[Plot[E^x+x^3/3+c1+x*c2,{x,-5,5},DisplayFunction->Identity],{c1,-10,10,5},{c2,-5,5,5}];Show[g1,DisplayFunction->$DisplayFunction];則輸出積分曲線的圖形(圖).圖例(教材例(教材例求微分方程組〈|dt+x+2y=et在初始條件x=1,y=0下的特解.xy=0t=0t=0輸入Clear[x,y,t];DSolve[{x'[t]+x[t]+2y[t]==Exp[t],y'[t]-x[t]-y[t]==0,x[0]==1,y[0]==0},{x[t],y[t]},t]例驗證1(5x230y+3y5)=c是微分方程y(x)=x2的通解.15y42輸入命令<<Graphics`PlotField`<<Graphics`ImplicitPlot`sol=(-5x^3-30y+3y^5)/15==C;g1=ImplicitPlot[sol/.Table[{C->n},{n,-3,3}],{x,-3,3}];g2=PlotVectorField[{1,x^2/(y^4-2)},{x,-3,3},{y,-3,3},Frame->True,ScaleFunction->(1&),ScaleFactor->,HeadLength->,PlotPoints->{20,25}];g=Show[g2,g1,Axes->None,Frame->True];Show[GraphicsArray[{g1,g2,g}]];則分別輸出積分曲線如圖(a),微分方程的方向場如圖(b).以及在同一坐標系中畫出積分曲線和方向場的圖形如下圖(c).圖從圖(c)中可以看出微分方程的積分曲線與方向場的箭頭方向吻合,且當x時,無論初始條件是什么,所有的解都趨向于一條直線方程.例(教材例求解微分方程dy2y=(x+1)5/2,并作出積分曲線.dxx+1輸入<<Graphics`PlotField`DSolve[y'[x]-2y[x]/(x+1)==(x+1)^(5/2),y[x],x]則輸出所給積分方程的解為下面在同一坐標系中作出這個微分方程的方向場和積分曲線(設輸入t=Table[2(1+x)^(7/2)/3+(1+x)^2c,{c,-1,1}];g1=Plot[Evaluate[t],{x,-1,1},PlotRange->{{-1,1},{-2,2}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];g2=PlotVectorField[{1,-2y/(x+1)+(x+1)^(5/2)},{x,,1},{y,-4,4},Frame->True,ScaleFunction->(1&),ScaleFactor->,HeadLength->,PlotPoints->{20,25},DisplayFunction->Identity];Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];則輸出積分曲線的圖形(圖).圖例求解微分方程(12xy)y=x2+y22,并作出其積分曲線.輸入命令<<Graphics`PlotField`DSolve[1-2*x*y[x]*y'[x]==x^2+(y[x])^2-2,y[x],x]則得到微分方程的解為t1=Table[(3+Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}];t2=Table[(3-Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}];gg1=Plot[Evaluate[t1],{x,-3,3},PlotRange->{{-3,3},{-3,3}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];gg2=Plot[Evaluate[t2],{x,-3,3},PlotRange->{{-3,3},{-3,3}},PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity];g1=ContourPlot[y-x^3/3-x*(-2+y^2),{x,-3,3},{y,-3,3},PlotRange->{-3,3},Contours->7,ContourShading->False,PlotPoints->50,DisplayFunction->Identity];g2=PlotVectorField[{1,(x^2+y^2-2)/(1-2*x*y)},{x,-3,3},{y,-3,3},Frame->True,ScaleFunction->(1&),ScaleFactor->,HeadLength->,PlotPoints->{20,25},DisplayFunction->Identity];Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];Show[gg1,gg2,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction];則輸出微分方程的向量場與積分曲線,并輸出等值線的圖.圖NDSolve似解問題:(1+xy)y+(1xy)y=0,y=1在區(qū)間[,4]上的近似解并作x=1.2.輸入fl=NDSolve[{(1+x*y[x])*y[x]+(1-x*y[x])*y'[x]==0,y[]==1},y,{x,,4}]則輸出為數(shù)值近似解(插值函數(shù))的形式:{{y->InterpolatingFunction[{{,4.}},<>]}}用Plot命令可以把它的圖形畫出來.不過還需要先使用強制求值命令Evalu-ate,輸入Plot[Evaluate[y[x]/.fl],{x,,4}]則輸出近似解的圖形(圖).圖如果要求區(qū)間[,4]內(nèi)某一點的函數(shù)的近似值,例如y,只要輸入x=1.8y[]/.fl則輸出所求結(jié)果例(教材例求范德波爾(VanderPel)方程在區(qū)間[0,20]上的近似解.輸入Clear[x,y];NDSolve[{y''[x]+(y[x]^2-1)*y'[x]+y[x]==0,y[0]==0,y'[0]==},y,{x,0,20}];Plot[Evaluate[y[x]/.%],{x,0,20}]可以觀察到近似解的圖形(圖).圖的數(shù)值解,并作出數(shù)值解的圖形.輸入命令<<Graphics`PlotField`sol=NDSolve[{x*y'[x]-x^2*y[x]*Sin[x]+1==0,y[1]==1},y[x],{x,1,4}];f[x_]=Evaluate[y[x]/.sol];g1=Plot[f[x],{x,1,4},PlotRange->All,DisplayFunction->Identity];g2=PlotVectorField[{1,(x^2*y*Sin[x]-1)/x},{x,1,4},{y,-2,9},Frame->True,ScaleFunction->(1&),ScaleFactor->,HeadLength->,PlotPoints->{20,25},DisplayFunction->Identity];g=Show[g1,g2,Axes->None,Frame->True];Show[GraphicsArray[{g1,g}],DisplayFunction->$DisplayFunction];則輸出所給微分方程的數(shù)值解及數(shù)值解的圖.圖例(教材例求出初值問題的數(shù)值解,并作出數(shù)值解的圖形.輸入NDSolve[{y''[x]+Sin[x]^2*y'[x]+y[x]==Cos[x]^2,y[0]==1,y'[0]==0},y[x],{x,0,10}]Plot[Evaluate[y[x]/.%],{x,0,10}];則輸出所求微分方程的數(shù)值解及數(shù)值解的圖形(圖).圖例(教材例洛倫茲(Lorenz)方程組是由三個一階微分方程組成的方程組.這三個方程看似簡單,也沒有包含復雜的函數(shù),但它的解卻很有趣和耐人尋味.試求解洛倫茲方程組并畫出解曲線的圖形.輸入Clear[eq,x,y,z]eq=Sequence[x'[t]==16*y[t]-16*x[t],y'[t]==-x[t]*z[t]-y[t]+45x[t],z'[t]==x[t]*y[t]-4z[t]];sol1=NDSolve[{eq,x[0]==12,y[0]==4,z[0]==0},{x[t],y[t],z[t]},{t,0,16},MaxSteps->10000];g1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol1],{t,0,16},PlotPoints->14400,Boxed->False,Axes->None];則輸出所求數(shù)值解的圖形(圖(a)).從圖中可以看出洛倫茲微分方程組具有一個奇異吸引子,這個吸引子緊緊地把解的圖形“吸”在一起.有趣的是,無論把解的曲線畫得多長,這些曲線也不相交.圖sol2=NDSolve[{eq,x[0]==6,y[0]==-10,z[0]==10},{x[t],y[t],z[t]},{t,0,24},MaxSteps->10000];g2=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol2],{t,0,24},PlotPoints->14400,Boxed->False,Axes->None];Show[GraphicsArray[{g1,g2}]];則輸出所求數(shù)值解的圖形(圖(b)).從圖中可以看出奇異吸引子又出現(xiàn)了,它把解“吸”在某個區(qū)域內(nèi),使得所有的解好象是有規(guī)則地依某種模式纏繞.實驗3拋射體的運動(續(xù))(綜合實驗)實驗目的通過微分方程建模和Mathematica軟件,在項目一實驗5的基礎上,進一步研究在考慮空氣阻力的情況下拋射體的運動.問題根據(jù)偵察,發(fā)現(xiàn)離我軍大炮陣地水平距離10km的前方有一敵軍的坦克群正以每小時50km向我軍陣地駛來,現(xiàn)欲發(fā)射炮彈摧毀敵軍坦克群.為在最短時間內(nèi)有效摧毀敵軍坦克,要求每門大炮都能進行精射擊,這樣問題就可簡化為單門大炮對移動坦克的精確射擊怎樣的發(fā)射角度可以最有效摧毀敵軍坦克.說明本節(jié)我們研究受到重力和空氣阻力約束的拋射體的射程.用(x(t),y(t))記拋計算出使拋射體的射程最大的發(fā)射角.假設t=0時拋射體(炮彈)在原點(0,0)以與水平線夾角為a,初始速度為v發(fā)射出去.它受到的空氣阻力為0F=_kv=_k(|dx,dy)|.r(dtdt)重力為gg在推導x(t)和y(t)所滿足的微分方程之前,補充一點說明:雖然我們將位置變量x(t),y(t)僅寫作t的函數(shù),但實際上位置變量還依賴于幾個其它的變量或參數(shù).特別ak、質(zhì)量m及重力加速度g等.為了推導x和y的方程,按照牛頓定律F=ma,并結(jié)合重力的公式和空氣阻力的公式,得到微分方程根據(jù)前面所述假設知,x(t),y(t)滿足下列初始條件x(0)=0,y(0)=0,x,(0)=vcosa,y,(0)=vsina.00先求解x(t),由方程,令v=x,,可將其化為一階微分方程易求出其通解由x(0)=0,得D=(|m)|vcosa,所以(k(k)x(t)=(|m)|vcosa(1_e_t)(k(k)類似地,可從方程解出y.令v=y,,方程化為一階微分方程,兩端除以m,得再在上述方程兩端乘以積分因子et.得即d(vet)=_get,dt兩端積分得vetgmetC.k所以vgmCet.k利用初始條件y(0)v(0)vsin確定其中的常數(shù)C后,積分v得到y(tǒng),再次利用初0始條件y(0)0確定任意常數(shù)后,則得到kkkk0下面我們利用公式與來描繪炮彈運行的典型圖形.假定炮彈發(fā)射的初速度為s,發(fā)射角為55。,輸入Clear[a,t,x,y,g,m,k]x[v_,a_,t_]:=(m/k)*v*Cos[aPi/180]*(1-Exp[-(k/m)*t])kmtmtg=;m=;k=;炮彈飛行的時間由炮彈落地時的條件y0所確定.輸入FindRoot[y[350,55,t]==0,{t,50}]則輸出炮彈飛行的時間x[350,55,]輸入ParametricPlot[{x[350,55,t],y[350,55,t]},{t,0,},PlotRange->{0,11000},AxesLabel->{x,y}]則輸出圖.圖實驗報告在上述假設下,進一步研究下列問題:(1)選擇一個初始速度和發(fā)射角,利用Mathematica畫出炮彈運行的典型軌跡.km的初速度為s,應選擇什么樣的發(fā)射角才能擊中坦克畫出炮彈運行的幾個軌跡圖,通過實驗數(shù)據(jù)和圖形來說明你的合理性.(3)假定坦克在大炮前方10km處靜止不動,探索降低或調(diào)高炮彈發(fā)射的初速度下,應如何選擇炮彈的發(fā)射角從上述討論中總結(jié)出最合理有效的發(fā)射速度和發(fā)射角.(4)在上題結(jié)論的基礎上,繼續(xù)探索,假定坦克在大炮前方10km處以每小時方向前進,此時應如何制定迅速摧毀敵軍坦克的方案化.實驗4蹦極跳運動(綜合實驗)實驗目的利用Mathematica軟件,通過微分方程建模,研究蹦極跳運動.問題在不考慮空氣阻力和考慮空氣阻力等多種情況下,研究蹦極跳運動中,蹦極者與蹦極繩設計之間的各種關(guān)系.說明蹦極繩相當于一根粗橡皮筋或有彈性的繩子.當受到張力使之超過其自然長度,繩子會產(chǎn)生一個線性回復力,即繩子會產(chǎn)生一個力使它恢復到自然長度,而這個力大小與它被拉伸的長度成正比.在一次完美的蹦極跳過程中,蹦極者爬上一座高橋或高的建筑物,把繩的一頭系在自己身上,另一頭系在一個固定物體如橋欄桿上,當他跳離橋時,激動人心的時刻就到來了.這里要分析的是蹦極者從跳出那一瞬間起他的運動規(guī)律.首先要建立坐標系.假設蹦極者的運動軌跡是垂直的,因此我們只要用一個坐標來確平面,時間tt=0,則y(t)表示t時刻蹦極者的位置.下面我們要求出y(t)的表達式.由牛頓第二定律,物體的質(zhì)量乘以加速度等于物體所受的力.我們假設蹦極者所受的力只有重力、空氣阻力和蹦極繩產(chǎn)生的回復力.當然,直到蹦極者降落的距離大于蹦極繩的自然長度時,蹦極繩才會產(chǎn)生回復力.為簡單起見,假設空氣阻力的大小與速度成正比,比例系數(shù)為1,蹦極繩回復力