# paramath

## * B1

Numerical and Graphical Solutions of Duffing Equation

by CO.H . TRAN  &   PHONG . T . NGO  -    University of Natural Sciences  , HCMC  Vietnam -
coth123@math.com   &  coth123@yahoo.com
Copyright  2006
Jan  06  2006
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
** Abstract  : Duffing quation is solved  by  Runge-Kutta  approximating method  .
** Subjects: Vibration Mechanics , The Differential equation  .
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

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

All rights reserved.  Copying or transmitting of this material without the permission of the authors is not allowed .                                       **********LỜI GIẢI SỐ VÀ ĐỒ THỊ PHƯƠNG TRÌNH DUFFING**********

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;

 >

 >

 > 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

Today, there have been 6 visitors (37 hits) on this page!
=> Do you also want a homepage for free? Then click here! <=