Search code examples
maple

How to obtain and work with steady-state part of DE solution?

I'm trying to obtain resonance curves of the system. System can be described as

F,m,k:=2,1,4: 
lambda:= beta/(2*m): 
omega:=sqrt(k/m):
de:=diff(x(t),t$2)+2*lambda*diff(x(t),t)+omega^2*x(t)=F*cos(gamma1*t):
cond:=x(0)=0, D(x)(0)=0:
sol := dsolve({cond, de});

Solving gives sum of terms, some of which "die out" with time (since these terms have exp(-...*t)) and some of which form steady-state solution (solution for t -> ∞). This solution will be in form xstst=f(gamma1)*sin(...). In order to obtain resonance curves, I need to plot f(gamma1) (for chosen constant betas, say, 2,1,0.5,0.25,etc.).

I've solved this "by hand" and found f := F/(sqrt((-gamma1^2+omega^2)^2+4*lambda^2*gamma1^2)). Plotting this for any chosen beta gives the result needed, for example, for beta:=0.5 the plot is enter image description here

I wonder if I can obtain these curves using maple functions only (without solving anything "by hand" at all).

[edited]enter image description here

Solution

  • It does not make sense to expect Maple to represent the result in terms of either sin(theta) or cos(theta), using some formulas for those terms that appear nowhere in the problem specification and are entirely introduced by you.

    The following is obtained using a radical (square root) in each of cond1 and cond2.

    restart;
de := diff(x(t),t,t)+2*lambda*diff(x(t),t)+omega^2*x(t)=F*cos(gamma1*t):
cond := x(0)=0, D(x)(0)=0:

sol := dsolve({cond, de}):

E,T := selectremove(hastype,rhs(sol),specfunc(anything,exp)):

lprint(T);

    F*(2*sin(gamma1*t)*gamma1*lambda-cos(gamma1*t)*gamma1^2
    +cos(gamma1*t)*omega^2)/(gamma1^4+4*gamma1^2*lambda^2-
    2*gamma1^2*omega^2+omega^4)

cond1 := cos(theta) = (-gamma1^2+omega^2)
            /((-gamma1^2+omega^2)^2+4*lambda^2*gamma1^2)^(1/2):

cond2 := sin(theta) = (-2*lambda*gamma1)
            /((-gamma1^2+omega^2)^2+4*lambda^2*gamma1^2)^(1/2):

frontend(algsubs, [numer(rhs(cond1))=lhs(cond1)*denom(rhs(cond1)),
                   numer(T)],
         [{`+`,`*`,`=`},{}]):

frontend(algsubs, [numer(rhs(cond2))=lhs(cond2)*denom(rhs(cond2)),
                   %],
         [{`+`,`*`,`=`},{}]):

ans := collect(combine(%, trig),cos)/denom(T):

lprint(ans);

    F*cos(gamma1*t+theta)/(gamma1^4+4*gamma1^2*lambda^2-
    2*gamma1^2*omega^2+omega^4)^(1/2)

subsans := eval(eval(ans,[lambda=beta/(2*m),omega=sqrt(k/m)]),
                [F=2,m=1,k=4]):

lprint(subsans);

    2*cos(gamma1*t+theta)
    /(beta^2*gamma1^2+gamma1^4-8*gamma1^2+16)^(1/2)