Как получить и работать с установившейся частью решения DE? - PullRequest
Новая
       10

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

0 голосов
/ 04 января 2019

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

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 график составляет enter image description here

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

[отредактировано] enter image description here

Ответы [ 2 ]

0 голосов
/ 31 января 2019

Нет смысла ожидать, что 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)
0 голосов
/ 17 января 2019

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

coeff( collect( rhs( sol ) , sin( gamma1 * t ) ) , sin( gamma1 * t ) , 1 )
...