paramath

* B1

Image 
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);
 

m*(diff(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 ;
 

y(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);
 

10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)-9000000*y(t)^3 = 100*cos(Pi*t)
10*(diff(diff(y(t), t), t))+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);
 

proc (x_rkf45) local res, data, vars, solnproc, outpoint, ndsol, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 14; if _EnvInF... 
proc (t) options operator, arrow; rhs(loigiai(t)[2]) end proc 
5 
[t = 0., y(t) = .85000000000000, diff(y(t), t) = 0.] 
[t = .20000000000000, y(t) = -0.192743776969386379e-4, diff(y(t), t) = -0.189311978181220356e-1]
[t = .20000000000000, y(t) = -0.192743776969386379e-4, diff(y(t), t) = -0.189311978181220356e-1]
 
[t = .40000000000000, y(t) = 0.309430997202991374e-5, diff(y(t), t) = -0.314943550511259446e-4]
[t = .40000000000000, y(t) = 0.309430997202991374e-5, diff(y(t), t) = -0.314943550511259446e-4]
 
[t = .60000000000000, y(t) = -0.308741697915723870e-5, diff(y(t), t) = -0.298814748639385770e-4]
[t = .60000000000000, y(t) = -0.308741697915723870e-5, diff(y(t), t) = -0.298814748639385770e-4]
 
[t = .80000000000000, y(t) = -0.808842606692218408e-5, diff(y(t), t) = -0.186448661419193254e-4]
[t = .80000000000000, y(t) = -0.808842606692218408e-5, diff(y(t), t) = -0.186448661419193254e-4]
 
[t = 1., y(t) = -0.100001048405655681e-4, diff(y(t), t) = 0.221276548367689596e-6]
[t = 1., y(t) = -0.100001048405655681e-4, diff(y(t), t) = 0.221276548367689596e-6]
 
Plot 
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;
 

5 
Plot 
Plot 
Plot 
Plot 
Plot 

>  
 

>  
 

> 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) ;
 

 Cac tham so : m =  
 Phuong trinh vi phan co dang : 
10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)+9000000*y(t)^3 = 100*cos(Pi*t)
10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)+9000000*y(t)^3 = 100*cos(Pi*t)
 
 Dieu kien dau :  
 y(0) =  
 y'(0) =  
 Loi giai so bang phuong phap RUNGE - KUTTA :  
[t = 0., y(t) = .50000000000000, diff(y(t), t) = 0.] 
[t = 1., y(t) = -0.100001922036141262e-4, diff(y(t), t) = 0.200259407474319750e-6]
[t = 1., y(t) = -0.100001922036141262e-4, diff(y(t), t) = 0.200259407474319750e-6]
 
[t = 2., y(t) = 0.100005096748307881e-4, diff(y(t), t) = 0.115874619398825230e-6]
[t = 2., y(t) = 0.100005096748307881e-4, diff(y(t), t) = 0.115874619398825230e-6]
 
[t = 3., y(t) = -0.999979035875673534e-5, diff(y(t), t) = 0.311729864059405209e-6]
[t = 3., y(t) = -0.999979035875673534e-5, diff(y(t), t) = 0.311729864059405209e-6]
 
[t = 4., y(t) = 0.999991013686392712e-5, diff(y(t), t) = -0.503447765508274115e-7]
[t = 4., y(t) = 0.999991013686392712e-5, diff(y(t), t) = -0.503447765508274115e-7]
 
[t = 5., y(t) = -0.100001052242845409e-4, diff(y(t), t) = -0.163032715969934755e-6]
[t = 5., y(t) = -0.100001052242845409e-4, diff(y(t), t) = -0.163032715969934755e-6]
 
Plot 

> mohinh(0.8);
 

Plot 

> mohinh(6);
 

Plot 

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

 Cac tham so : m =  
 Phuong trinh vi phan co dang : 
10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)+9000000*y(t)^3 = 100*cos(Pi*t)
10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)+9000000*y(t)^3 = 100*cos(Pi*t)
 
 Dieu kien dau :  
 y(0) =  
 y'(0) =  
 Loi giai so bang phuong phap RUNGE - KUTTA :  
[t = 0., y(t) = .15000000000000, diff(y(t), t) = 0.] 
[t = 1., y(t) = -0.100002100965684427e-4, diff(y(t), t) = 0.174445319969807240e-6]
[t = 1., y(t) = -0.100002100965684427e-4, diff(y(t), t) = 0.174445319969807240e-6]
 
[t = 2., y(t) = 0.100003219245372214e-4, diff(y(t), t) = -0.758163623338622300e-7]
[t = 2., y(t) = 0.100003219245372214e-4, diff(y(t), t) = -0.758163623338622300e-7]
 
[t = 3., y(t) = -0.100002374849772969e-4, diff(y(t), t) = -0.190236480746275206e-6]
[t = 3., y(t) = -0.100002374849772969e-4, diff(y(t), t) = -0.190236480746275206e-6]
 
[t = 4., y(t) = 0.100004144928036077e-4, diff(y(t), t) = 0.770508220812301378e-7]
[t = 4., y(t) = 0.100004144928036077e-4, diff(y(t), t) = 0.770508220812301378e-7]
 
[t = 5., y(t) = -0.100003957106223904e-4, diff(y(t), t) = -0.890413129754806662e-7]
[t = 5., y(t) = -0.100003957106223904e-4, diff(y(t), t) = -0.890413129754806662e-7]
 
[t = 6., y(t) = 0.100002617106955826e-4, diff(y(t), t) = -0.153364006735818149e-6]
[t = 6., y(t) = 0.100002617106955826e-4, diff(y(t), t) = -0.153364006735818149e-6]
 
[t = 7., y(t) = -0.100000618081845140e-4, diff(y(t), t) = 0.264609326752723164e-6]
[t = 7., y(t) = -0.100000618081845140e-4, diff(y(t), t) = 0.264609326752723164e-6]
 
[t = 8., y(t) = 0.99999993278359298e-5, diff(y(t), t) = -0.130065877267496021e-6]
[t = 8., y(t) = 0.99999993278359298e-5, diff(y(t), t) = -0.130065877267496021e-6]
 
[t = 9., y(t) = -0.99999913174129808e-5, diff(y(t), t) = 0.824227836181808212e-7]
[t = 9., y(t) = -0.99999913174129808e-5, diff(y(t), t) = 0.824227836181808212e-7]
 
[t = 10., y(t) = 0.999991202646871412e-5, diff(y(t), t) = -0.171616608259521394e-6]
[t = 10., y(t) = 0.999991202646871412e-5, diff(y(t), t) = -0.171616608259521394e-6]
 
Plot 

> mohinh(5);
 

Plot 

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);
 

 Cac tham so : m =  
 Phuong trinh vi phan co dang : 
10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)+9000000*y(t)^3 = 100*cos(Pi*t)
10*(diff(diff(y(t), t), t))+1000*(diff(y(t), t))+10000000*y(t)+9000000*y(t)^3 = 100*cos(Pi*t)
 
 Dieu kien dau :  
 y(0) =  
 y'(0) =  
 Loi giai so bang phuong phap RUNGE - KUTTA :  
[t = 0., y(t) = 0., diff(y(t), t) = .10000000000000] 
[t = 1., y(t) = -0.100003031557671699e-4, diff(y(t), t) = 0.856054697571733248e-7]
[t = 1., y(t) = -0.100003031557671699e-4, diff(y(t), t) = 0.856054697571733248e-7]
 
[t = 2., y(t) = 0.100001186730557259e-4, diff(y(t), t) = 0.138567705758366402e-6]
[t = 2., y(t) = 0.100001186730557259e-4, diff(y(t), t) = 0.138567705758366402e-6]
 
[t = 3., y(t) = -0.100001881829770924e-4, diff(y(t), t) = 0.143173631414504646e-6]
[t = 3., y(t) = -0.100001881829770924e-4, diff(y(t), t) = 0.143173631414504646e-6]
 
[t = 4., y(t) = 0.999992373954906582e-5, diff(y(t), t) = 0.660995478059635895e-7]
[t = 4., y(t) = 0.999992373954906582e-5, diff(y(t), t) = 0.660995478059635895e-7]
 
[t = 5., y(t) = -0.100003314492074806e-4, diff(y(t), t) = 0.314077500376234288e-7]
[t = 5., y(t) = -0.100003314492074806e-4, diff(y(t), t) = 0.314077500376234288e-7]
 
[t = 6., y(t) = 0.999992793437681750e-5, diff(y(t), t) = 0.153023326407566206e-6]
[t = 6., y(t) = 0.999992793437681750e-5, diff(y(t), t) = 0.153023326407566206e-6]
 
[t = 7., y(t) = -0.100001940930439043e-4, diff(y(t), t) = -0.173254913410583658e-6]
[t = 7., y(t) = -0.100001940930439043e-4, diff(y(t), t) = -0.173254913410583658e-6]
 
[t = 8., y(t) = 0.100005562067931500e-4, diff(y(t), t) = -0.219553498185131744e-6]
[t = 8., y(t) = 0.100005562067931500e-4, diff(y(t), t) = -0.219553498185131744e-6]
 
[t = 9., y(t) = -0.100009141784075014e-4, diff(y(t), t) = -0.327149509524763626e-6]
[t = 9., y(t) = -0.100009141784075014e-4, diff(y(t), t) = -0.327149509524763626e-6]
 
[t = 10., y(t) = 0.100003112333562182e-4, diff(y(t), t) = -0.682392457566916905e-7]
[t = 10., y(t) = 0.100003112333562182e-4, diff(y(t), t) = -0.682392457566916905e-7]
 
Plot 

> mohinh(5);
 

Plot 

>  
 

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

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.  

Image 
B1.COHONGTRAN-DUFFEQN








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