Как получить и работать с установившейся частью решения ДУ?

Я пытаюсь получить кривые резонанса системы. Система может быть описана как

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

Решение дает сумму слагаемых, часть которых со временем "вымирает" (поскольку эти слагаемые имеют exp(-...*t)), а часть образует стационарное решение (решение для t -> ∞). Это решение будет иметь форму xstst=f(gamma1)*sin(...). Чтобы получить кривые резонанса, мне нужно построить f(gamma1) (для выбранной константы betas, скажем, 2,1,0.5,0.25,и т. д.).

Я решил это "вручную" и нашел f := F/(sqrt((-gamma1^2+omega^2)^2+4*lambda^2*gamma1^2)). Построение этого графика для любой выбранной бета-версии дает необходимый результат, например, для beta:=0.5 график введите здесь описание изображения

Интересно, смогу ли я получить эти кривые, используя только функции Maple (вообще ничего не решая «вручную»).

[отредактировано]введите здесь описание изображения


maple
person Kelly Shepphard    schedule 03.01.2019    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Не имеет смысла ожидать, что Maple представит результат в терминах sin(theta) или cos(theta), используя некоторые формулы для этих терминов, которые нигде не встречаются в описании задачи и полностью введены вами.

Следующее получается с использованием радикала (квадратный корень) в каждом из cond1 и 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)
person acer    schedule 31.01.2019

arrow_upward
1
arrow_downward

Я не уверен, что полностью понимаю ваш вопрос, но когда я запускаю ваш код, я получаю некоторые члены с экспонентами, некоторые с синусом и некоторые с косинусами. Вы можете получить коэффициент синуса с помощью

coeff( collect( rhs( sol ) , sin( gamma1 * t ) ) , sin( gamma1 * t ) , 1 )
person Ian Thompson    schedule 16.01.2019
comment
Спасибо за ответ, но мне нужно найти coeff только для последнего слагаемого (которое не содержит exp(-t*...)) и терм будет как f(gamma1)*sin(gamma1*t+theta ), где тета также является функцией лямбда, омега и гамма1. Я пытался собрать это, но кажется, что это не сработает (собирание дает тот же 3-й член, что и раньше). По крайней мере, теперь я знаю, что это невозможно. +1 за помощь, спасибо. Эта идея оказалась очень полезной для понимания. - person Kelly Shepphard; 17.01.2019
comment
@KellyShepphard --- Я подозреваю, что это возможно, но я все еще не уверен на 100%, что именно вы хотите. (В вашем коде нет теты, и я не знаю, откуда берется sqrt в работе, которую вы проделали вручную.) Возможно, я смогу дать лучший ответ, если вы загрузите скриншот из Maple и выделите нужные термины. изолировать. - person Ian Thompson; 17.01.2019
comment
Я добавил скриншот. Извините за отсутствие информации, я думал, что Maple автоматически решает ее по мере необходимости. - person Kelly Shepphard; 17.01.2019