این برای درس موج در زمستان سال ۹۶ آماده شده است. در این بخش به موضوع نوسانات واداشته می پردازیم.

پیش از شروع

SetDirectory[NotebookDirectory[]];
PlotDir=Directory[]<>"/Plots_Data_SOMQ/";(*The path to store the plots.*)

If[DirectoryQ[PlotDir],Print["Is already there"],CreateDirectory[PlotDir]]

Exporter[name_]:=Export[PlotDir<>name,%];

fontS=24;imgS=Large;Themes`AddThemeRules["SadeghNicePlots",(*Default*)PlotStyle->Thread@Directive[{Red,Blue,Purple,Orange,Hue[0.6]},Thickness[.008]],LabelStyle->fontS,LabelStyle->Directive[Bold, FontSize-> fontS],ImageSize->imgS,AxesStyle->{{Bold,Thick,FontSize-> 20},{Bold,Thick,FontSize-> 20}},TicksStyle->{{Gray,Thick},{Gray,Thick}},Background->White,PlotRange->All,PlotMarkers->{Automatic,10}]

$PlotTheme=SadeghNicePlots

$PlotTheme="Default"

تعاریف مقدماتی

m=1;k=1;Γ=.1;fd=.5;ωd=2;ω0=(Sqrt[(k/m)]);x0=1;v0=1;t0=100*2π/ωd;

Sqrt[2 ω0]

Ftotal[t_,x_,v_, m_,k_,Γ_,f0_,ωd_]:=- Γ m  v-k x+f0 Cos[ωd t];Ftotal[t,x0,v0,M121,kspring,γ,f0,Ωd]

Ftotal[t,x0 Cos[Ωd t] ,-Ωd v0 Sin[Ωd t], M,Ks,Γs,f0,Ωd]

Plot[Ftotal[t,x0 Cos[ωd t] ,-ωd v0 Sin[ωd t], m,k,Γ,fd,ωd], {t,0,3*2*π/ωd}]

Ftotal[t,x[t],x'[t],m,k,Γ,fd,ωd]

{time,sol}=Timing[DSolve[{m x''[t]==Ftotal[t,x[t],x'[t],m,k,Γ,fd,ωd],x[0]==.1x0,x'[0]==0},x[t],t]]

X[t_] = (x[t]/.sol)[[1]];X[1]//Chop

positionList  = Table[{t,Chop[X[t]]},{t,0,2π/ωd,2π/ωd/100}];ListPlot[positionList]

t0=100*2π/ωd;positionList  = Table[{t,Chop[X[t]]},{t,t0,t0+2π/ωd,2π/ωd/100}];ListPlot[positionList]

(*positionList  = Table[{t,Chop[X[t]]},{t,0,t0+2π/ωd,π/ωd/7}];ListPlot[positionList,Joined->True]*)Plot[X[t],{t,0,t0},PlotPoints->1000,WorkingPrecision->20]

Clear@sol

{time,sol[t_]}=Timing[NDSolve[{m x''[t]==Ftotal[t,x[t],x'[t],m,k,Γ,fd,ωd],x[0]==0.1,x'[0]==0},x[t],{t,0,t0+2π/ωd}]]

Pos[t_]:=(x[t]/.sol[t])[[1]]

Vel[t_]:=Evaluate[D[Pos[t],t]]

Pos[1]

Vel[1]

Plist= Table[{Pos[t],Vel[t]},{t,0,10π/ωd,10π/ωd/100}];ListPlot[Plist,Joined->True]

Plot[Evaluate[x[t]/.sol],{t,t0,t0+2π/ωd},PlotRange->All]

Plot[Evaluate[x[t]/.sol],{t,0,t0/4+2π/ωd},PlotRange->All]

جواب های پایا

Clear[ω0,ωd,Γ,sol]

sol[ωd_]= Solve[A(((ω0)^2)-((ωd)^2)-I ωd Γ) == f0/M,A];A[ωd_]:=(A/.sol[ωd])[[1]];A[Ω]

With[{ω0=.8,Γ=.1,f0=1,M=m},Plot[{Re[(f0/(M (-I Γ Ω-(Ω^2)+((ω0)^2))))],Im[(f0/(M (-I Γ Ω-(Ω^2)+((ω0)^2))))]},{Ω,0,2},AxesLabel->{"ωd","Amp"},PlotLegends->{"A","B"}]]

Manipulate[Plot[{Re[(f0/(M (-I Γ Ω-(Ω^2)+((ω0)^2))))],Im[(f0/(M (-I Γ Ω-(Ω^2)+((ω0)^2))))]},{Ω,0,2},AxesLabel->{"ωd","Amp"},PlotLegends->{"A","B"}],{ω0,.1,2},{Γ,.1,1},{{f0,1},0,2},{{M,1},.1,2}]

گاهی اوقات نرم افزار توانایی درک برخی از خصوصیات را ندارد. در مثال زیر می بینید که نمی تواند قسمت موهومی تابع را ندارد. در حالی که ما می توانی به سادگی این را حساب کنیم. .

FullSimplify[Im[(f0/(M (-I Γ Ω-(Ω^2)+((ω0)^2))))],f0>0&&M>0&&Ω>0&&ω0>0&&Γ>0]

به همین دلیل گاهی بهتر است خودمان توابع را مستقیم تعریف کنیم.

A[Ω_,ω0_,Γ_,f0_,M_]:=(f0/(M ))((-(Ω^2)+((ω0)^2))/((((( Γ Ω))^2)+(((-(Ω^2)+((ω0)^2)))^2))));
B[Ω_,ω0_,Γ_,f0_,M_]:=(f0/(M ))((Γ Ω)/((((( Γ Ω))^2)+(((-(Ω^2)+((ω0)^2)))^2))));

Clear[m,Γ,fd,k,ω0,ωd]

X[t_,Ω_,ω0_,Γ_,f0_,M_]:=(A [Ω,ω0,Γ,f0,M]Cos[Ω t]+B [Ω,ω0,Γ,f0,M]Sin[Ω t]);
V[t_,Ω_,ω0_,Γ_,f0_,M_]:=Ω(-A [Ω,ω0,Γ,f0,M]Sin[Ω t]+B [Ω,ω0,Γ,f0,M]Cos[Ω t]);

می توان دید که مشتق همان تابع سرعت را می دهد.

D[(A Cos[Ω t]+B Sin[Ω t]),t]

تاخر فاز

در زیر می توانید مقایسه نوسانات نیروی خارجی را در کنار مکان سیستم ببینید. این تاخر فاز پاسخ را نسبت به نیروی خارجی نشان میدهد.

Manipulate[Plot[{f0 Cos[Ω t],A[Ω,ω0,Γ,f0,M]Cos[Ω t]+B[Ω,ω0,Γ,f0,M]Sin[Ω t]},{t,0,20},AxesLabel->{"Time","Amp"},PlotLegends->{"F","X"}],{{Ω,1},.1,2},{ω0,.1,2},{Γ,.1,1},{{f0,1},0,2},{{M,1},.1,2}]

نمودار زیر نیز اختلاف فاز بین نوسانات نیرو و مکان نوسناگر را نشان می دهد.

Manipulate[Plot[Arg[I B[Ω,ω0,Γ,f0,M]+ A[Ω,ω0,Γ,f0,M]],{Ω,0,2},AxesLabel->{"ωd","Phase Lag"},PlotRange->{All,{0,1.2π}},Ticks->{{0,ω0,2ω0},{0,π/2,π}}],{{ω0,1},.1,2},{Γ,.1,1},{{f0,1},0.1,2},{{M,1},.1,2}]

انرژی و توان

P[t_,Ω_,ω0_,Γ_,f0_,M_]:=f0(-A [Ω,ω0,Γ,f0,M]Cos[Ω t]Sin[Ω t]+B [Ω,ω0,Γ,f0,M]((Cos[Ω t])^2));

Manipulate[Plot[{P[t,Ω,ω0,Γ,f0,M],-f0 A [Ω,ω0,Γ,f0,M]Cos[Ω t]Sin[Ω t],f0 B [Ω,ω0,Γ,f0,M]((Cos[Ω t])^2)},{t,0,20},AxesLabel->{"Time","Power"},PlotLegends->{"Total Power","Part A","Part B"}],{{Ω,1},.1,2},{ω0,.1,2},{Γ,.1,1},{{f0,1},0,2},{{M,1},.1,2}]

دو قسمت نوسانی در توان وجود دارد. یکی از این بخش ها میانگینش صفر است و دیگری غیر صفر. انتگرال های زیر این را نشان می دهد.

Integrate[f0(-A Cos[Ω t]Sin[Ω t]),{t,0,2π/Ω}]

FullSimplify[Integrate[f0 Ω(B ((Cos[Ω t])^2)),{t,0,2π/Ω}],f0>0&&M>0&&Ω>0&&ω0>0&&Γ>0]

Integrate[f0 Ω(B ((Cos[Ω t])^2)),t]

Pavg[Ω_,ω0_,Γ_,f0_,M_]:=((Ω B[Ω,ω0,Γ,f0,M] f0 )/2)

Manipulate[Plot[Pavg[Ω,ω0,Γ,f0,M],{Ω,0,2},AxesLabel->{"ωd","Amp"}(*,PlotLegends->{"A","B"}*)],{{ω0,1},.1,2},{Γ,.1,1},{{f0,1},0,2},{{M,1},.1,2}]

Manipulate[Plot[Pavg[Ω,ω0,Γ,f0,M]+Pavg[Ω,ω1,Γ,f0,M],{Ω,0,2},AxesLabel->{"ωd","Amp"}(*,PlotLegends->{"A","B"}*)],{{ω0,1},.1,2},{Γ,.1,1},{{f0,1},0,2},{{ω1,1},.1,2},{{M,1},.1,2}]

Manipulate[Plot[Pavg[Ω,ω0,Γ,f0,M]+Pavg[Ω,ω1,Γ,f0,M]+Pavg[Ω,ω2,Γ,f0,M],{Ω,0,2},AxesLabel->{"ωd","Amp"}(*,PlotLegends->{"A","B"}*)],{{ω0,1},.1,2},{{ω1,1},.1,2},{Γ,.1,1},{{f0,1},0,2},{{ω2,1},.1,2},{{M,1},.1,2}]