paramath

* B1

 

                                                                                                                                                                                               

NOTE: This worksheet demonstrates Maple's  capabilities in the design of a dynamic system and finding  the numerical solution of the Duffing equation .                                     

                                                         TRAN HONG CO  &  NGO THANH PHONG  - Dai hoc Khoa hoc tu nhien - tp HCM  Vietnam 

                                                                                       coth123@math.com  &  coth123@yahoo.com   

A .  Xac dinh he thong .  [ System Definition ] 

>

restart:

 

>	with(plots):

 

>	interface(warnlevel=0):

 

B. Mo hinh dao dong phi tuyen .  [ Non-linear Vibration Model ] 

Khao sat mau vat the  dang giai tich co hinh mo phong nhu tren  [ Consider an analytical model which has the simulation figure above ] 

Phuong trinh vi phan chuyen dong : [ Differential equation of vibration ] 

>	ptcd:= m*diff(y(t),t,t) + b*diff(y(t),t) + c1*y(t) + c3*y(t)^3  =  F(t);

 

 

Xac dinh cac dieu kien dau  . [ Define initial conditions ] 

>	dkdau:= y(0)=0.85,D(y)(0)=0 ;

 

 

Thay luc F phi tuyen va cac gia tri cua tham so m , b , c1 , c3 .  [ Substitute non-linearity force  F and parameter values  m , b , c1 , c3 ] . 

>	ptcdcuthe:= 10*diff(y(t),`$`(t,2))+1000*diff(y(t),t)+10000000*y(t)-9000000*y(t)^3 = 100*cos(Pi*t);

 

 

Loi giai so cua phuong trinh vi phan phi tuyen  . [ Numerical solution of the non-linear differential equation ] 

>	loigiai:=dsolve({ptcdcuthe,dkdau}, y(t), type=numeric, method=rkf45);with(DEtools):with(plots):u:=t->rhs(loigiai(t)[2]);n:= 5;for i from 0 to n do    loigiaiso[i]:=normal(loigiai(i/n)); od;plot(u,-1..5,axes=box);

 

 

 

 

 

 

 

 

 

 

 

Do thi dao dong voi  t  =  0  (giay)  den  t  =  1  (giay)  [ Graph of vibration  from  t  =  0 (s)   to   t  =  1 (s) 

>	m:=5;for j from 1 to m do plot(u,j/(m)..0,axes=box); od;

 

 

 

 

 

 

 

>

 

C . Chuong trinh tinh toan mo hinh dao dong phi tuyen  .  [ Calculation procedures for the non-linear vibration model ] 

Bang hoat trinh Maple  nay gom 2 thu tuc . Ví du ve cach su dung trong cac bai toan thuc te , hay xem phan tiep sau . 

* daodong(khoi luong M , hang so can nhot b, he so lo xo c1, he so lo xo c3, luc F , dieu kien dau  , hoanh do) 

This Maple worksheet contains 2 procedures .  For examples of applying them to real problems, see the following action  .  

* daodong(mass M ,viscous damping constant b,spring constant c1, spring constant c3, force F , initial conditions ,absciss) 

 

>

 

>

daodong:=proc(Mf,bf,c1f,c3f,Ff,Nf,dkd,hd)                                                                                            global tungdo,y1,N;                                                                                                                   local m,b,c1,c3,k,F,omega0,g,tungdo0,daohamtungdo0,eq1,G;  with(DEtools): tungdo0:=y(0); daohamtungdo0:=D(y)(0);N:=Nf;m:=Mf;b:=bf;c1:=c1f;c3:=c3f;;F:=Ff;;; eq1:=m*diff(y(t),t,t) + b*diff(y(t),t) + c1*y(t) + c3*y(t)^3  =  F;print(" Cac tham so : m = ",m,"b = ",b,"c = ",c,"F = ",F,"N = ",N);print(" Phuong trinh vi phan co dang :");print(eq1);;print(" Dieu kien dau : ");print(" y(0) = ",y[0],"=",dkd[1]); print(" y'(0) = ",daohamtungdo0,"=",[2]); ;;G:=dsolve({eq1,tungdo0=dkd[1],daohamtungdo0=dkd[2]},{y(t)},type=numeric,method=rkf45); tungdo:=t-> rhs(G(t)[2]):print(" Loi giai so bang phuong phap RUNGE - KUTTA : ");for k from 0 to N do print(G(k)) od:;y1:=t->tungdo(t) : :;with(plottools):with(plots):plot([tungdo],hd+1/N..0,-0.0001..0.0001,color=[red],axes=frame,thickness=[1],title=`tung do mau do`);end: print("---------------------------------------------------------------------------------------------------------------------------------");                                                                                                                   mohinh:=proc(M) local Gi,i,j,omega,xt,yt,l,A,li,day,dia,thanh1a,thanh1b,thanh2,thanh3a,thanh3b,thanh4a,thanh4b,vatP,vatQ,loxo,cannhot,dem,nenlon,nennho,hinh,vanh1,vanh2,vanh3,vanh4,vanh5,vanh6,vanh7,tamgiac,vienbilon,vienbinho,bigiua,vong2a,vong2b,vong2c,vong2d,vong2e,tayquay,vo1a,vo1b,vo1c,lucF,lucFa,lucFb;with(plottools):Gi:=10*M;for i from 0 to Gi do thanh1a:=plot([[1,-0.5-y1(i/10)],[1,y1(i/10)-3]],style=line,thickness=3,color=black,view=[-10..10,-10..10]):thanh1b:=plot([[1,-5.5],[1,-7]],style=line,thickness=3,color=black,view=[-10..10,-10..10]):thanh2:=plot([[0,2.5-y1(i/10)],[0,5-y1(i/10)]],thickness=3,style=line,color=blue,view=[-10..10,-10..10]):thanh3a:=plot([[-1,-0.5-y1(i/10)],[-1,y1(i/10)-3.8]],style=line,thickness=3,color=grey,view=[-10..10,-10..10]):thanh3b:=plot([[-1,-4.75],[-1,-7]],style=line,thickness=3,color=grey,view=[-10..10,-10..10]): thanh4a:=plot([[-0.35,3.5-y1(i/10)],[-0.35,1.5-y1(i/10)]],thickness=7,style=line,color=black,view=[-10..10,-10..10]):thanh4b:=plot([[0.35,3.5-y1(i/10)],[0.35,1.5-y1(i/10)]],thickness=7,style=line,color=black,view=[-10..10,-10..10]): vienbilon:=plottools[circle]([0,6],1,1,color=black,view=[-1.5..1.5,-1.5..1.5]):vienbinho:=plottools[circle]([0,6],0.5,0.2,color=black,thickness=3,view=[-1.5..1.5,-1.5..1.5]): tayquay:=plot([[0,6],[-9+3*y1(i/10),+6+3*y1(i/10)]],thickness=5,style=line,color=black,view=[-10..10,-10..10]): vatP:=plottools[rectangle]([-3,-y1(i/10)+1.7],[3,-y1(i/10)-1.7],color=red,view=[-10..10,-10..10]):loxo:=plottools[rectangle]([-1.5,-3*y1(i/10)-4.5],[-0.5,-3*y1(i/10)-2.5],style=point,color=white,view=[-10..10,-10..10]): vatQ:=plottools[rectangle]([-0.25,5.5-y1(i/10)],[0.25,2-y1(i/10)],color=yellow,view=[-10..10,-10..10]): ;cannhot:=plottools[rectangle]([0.5,-5.5],[1.5,-2.5],color=yellow,view=[-10..10,-10..10]): dem:=plottools[rectangle]([0.5,-2*y1(i/10)-2.5],[1.5,-2*y1(i/10)-2.85],color=white,view=[-10..10,-10..10]):vanh1:=plottools[rectangle]([-1.5,-3*y1(i/10)-4.4],[-0.5,-3*y1(i/10)-4.3],color=black,view=[-10..10,-10..10]):vanh2:=plottools[rectangle]([-1.5,-3*y1(i/5)-3.7],[-0.5,-3*y1(i/5)-3.8],color=black,view=[-10..10,-10..10]):vanh3:=plottools[rectangle]([-1.5,-3*y1(i/5)-3.5],[-0.5,-3*y1(i/5)-3.4],color=black,view=[-10..10,-10..10]):vanh4:=plottools[rectangle]([-1.5,-3*y1(i/5)-2.4],[-0.5,-3*y1(i/5)-2.3],color=black,view=[-10..10,-10..10]):vanh5:=plottools[rectangle]([-1.5,-3*y1(i/5)-2.7],[-0.5,-3*y1(i/5)-2.8],color=black,view=[-10..10,-10..10]):vanh6:=plottools[rectangle]([-1.5,-3*y1(i/5)-4.5],[-0.5,-3*y1(i/5)-4.6],color=black,view=[-10..10,-10..10]):vanh7:=plottools[rectangle]([-1.5,-3*y1(i/5)-4.2],[-0.5,-3*y1(i/5)-4.1],color=black,view=[-10..10,-10..10]):nenlon:=plottools[rectangle]([-10,-8.5],[10,-10],color=green,view=[-10..10,-10..10]):;nennho:=plottools[rectangle]([-5,-6.95],[5,-8.5],color=blue,view=[-10..10,-10..10]): bigiua:=plottools[circle]([0,0-y1(i/10)],0.4,-0.4,color=yellow,thickness=3,view=[-1.5..1.5,-1.5..1.5]): vong2a:=plot([[-0.8,2.6-y1(i/10)],[0.8,2.6-y1(i/10)]],style=line,thickness=3,color=blue,view=[-10..10,-10..10]): vong2b:=plot([[-0.8,3.7-y1(i/10)],[0.8,3.7-y1(i/10)]],style=line,thickness=3,color=blue,view=[-10..10,-10..10]): vong2c:=plot([[-0.8,4.8-y1(i/10)],[0.8,4.8-y1(i/10)]],style=line,thickness=3,color=blue,view=[-10..10,-10..10]): vong2d:=plot([[-0.8,2.9-y1(i/10)],[0.8,2.9-y1(i/10)]],style=line,thickness=3,color=blue,view=[-10..10,-10..10]): vong2e:=plot([[-0.8,4.2-y1(i/10)],[0.8,4.2-y1(i/10)]],style=line,thickness=3,color=blue,view=[-10..10,-10..10]): :lucF:=plot([[0,8.5+y1(i/10)],[0,12-y1(i/10)]],thickness=5,style=line,color=red,view=[-10..10,-10..10]): :lucFa:=plot([[-0.35,10+y1(i/10)],[0,9-y1(i/10)]],thickness=5,style=line,color=red,view=[-10..10,-10..10]): :lucFb:=plot([[0,9-y1(i/10)],[0.35,10+y1(i/10)]],thickness=5,style=line,color=red,view=[-10..10,-10..10]): vo1a:=plot([[-10,-10],[-10,15]],style=line,thickness=3,color=black,view=[-10..10,-10..10]):vo1b:=plot([[10,-10],[10,15]],style=line,thickness=3,color=black,view=[-10..10,-10..10]):;vo1c:=plot([[-10,15],[10,15]],style=line,thickness=3,color=black,view=[-10..10,-10..15]):tamgiac:= polygon([[-0.35,3.5-y1(i/5)],[0,2.5-y1(i/5)],[0.35,3.5-y1(i/5)]], color=blue, thickness=3):l:=6.2;: A:=0.1; omega:=180;li:=5:;;;;;;;;for j from 0 by 1 to i do xt := l*sin(Pi*(j/li)): yt :=l*cos(Pi*(j/li)): day[j]:=curve([[0,-A*cos(omega*(i/li))], [xt,yt]],thickness=2,color=black):dia[j]:=disk([xt,yt], 0.6, color=black):end do:          hinh[i]:=plots[display]({thanh1a,thanh1b,thanh2,thanh3a,thanh3b,vatP,vatQ,cannhot,dem,nenlon,nennho,vanh1,vanh2,vanh3,vanh4,vanh5,vanh6,vanh7,thanh4a,thanh4b,tayquay,vienbilon,vienbinho,bigiua,vo1a,vo1b,vo1c,vong2a,vong2b,vong2c,vong2d,vong2e,lucF,lucFa,lucFb,dia[i],day[i]}):od:plots[display]([seq(hinh[i],i=0..Gi)],insequence=true,axes=box,scaling=constrained,view=[-1.2*l-10..1.2*l+10, -1.2*l-10..1.2*l+10],title=`lucF-do,M-do,b-vang,c1-xam,c3-xanh,TRANHONGCO`):end:

 

 

>	daodong(10,1000,10000000,9000000,100*cos(Pi*t),5,[0.5,0],5.8) ;

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

>	mohinh(0.8);

 

 

>	mohinh(6);

 

 

>	daodong(10,1000,10000000,9000000,100*cos(Pi*t),10,[0.15,0],1.8);

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

>	mohinh(5);

 

 

Gia su ta co dieu kien dau khong tuong thich  , dieu gi se xay ra ?   [ Suppose we have the incompatible initial conditions , what then  ? ] 

Ta co the thay ket qua bat thuong tu do thi va mo phong dao dong . [ We can see the extraordinary result from the  graph and simulation of vibration model ] 

>	daodong(10,1000,10000000,9000000,100*cos(Pi*t),10,[0,0.1],4.8);

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

>	mohinh(5);

 

 

>

 

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 

Disclaimer: While every effort has been made to validate the solutions in this worksheet, the authors  are  not responsible for any errors contained and are not liable for any damages resulting from the use of this material. 

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ 

Legal Notice: The copyright for this application is owned by Maplesoft. The application is intended to demonstrate the use of Maple to solve a particular problem. It has been made available for product evaluation purposes only and may not be used in any other context without the express permission of Maplesoft.  

 

B1.COHONGTRAN-DUFFEQN