# How to find and plot both steady periodic solution and actual solution?

Hi! I haven't been able to figure out how to do this in Sage, Let's say I have the differential equation with initial values

x''+4x'+5x=15cos(4t), x(0)=0, x'(0)=0


And want to find both the steady periodic solution as well as the full solution

When I put in

f=function('f')(x)
de=diff(f,x,2)+4*diff(f,x)+5*f==15*cos(4*x)
z=desolve_laplace(de,dvar=f,ics=[0,0,0])
show(expand(z))


I get

165/377*cos(x)*e^(-2*x) - 630/377*e^(-2*x)*sin(x) - 165/377*cos(4*x) + 240/377*sin(4*x)


Xsp(t) = 0.773cos(4t-2.173)
X(t) = 0.773cos(4t-2.173)+1.727e^(-2t) *cos(t-4.969)


And I have no idea how to graph this, but I'm pretty sure that's fairly easy. I know there's the manual way of doing everything, but that takes a lot of time and I think Sage can do it better and faster

I appreciate all help!

edit retag close merge delete

Sort by » oldest newest most voted

[ Note : you should try to avoit fixing the arguments a=of a symbolic function... ]

sage: f = function('f')
sage: de = diff(f(x),x,2) + 4*diff(f(x),x) + 5*f(x) == 15*cos(4*x)
sage: PM = desolve(de, f(x), ivar=x, ics=(0,0,0))
sage: PM
15/377*(11*cos(x) - 42*sin(x))*e^(-2*x) - 165/377*cos(4*x) + 240/377*sin(4*x)


The expressions of the shape $a\cos(u)+ b\sin(u)$ can indeed be expressed as $c\cos(u+d)$, but that is irrelevant to your question. More relevant is to consider your function as a polynomial of the dampening function $e^{-2x}$

sage: PM.subs(e^(-2*x)==E2X).coefficients(E2X)
[[-165/377*cos(4*x) + 240/377*sin(4*x), 0],
[165/377*cos(x) - 630/377*sin(x), 1]]


This decomposition allows to decompose your function as the sum of a "steady" and a "dampened" functions :

sage: fs, fd = ((u*E2X^u).subs(E2X==e^(-2*x)).function(x) for u in PM.subs(e^(-2*x)==E2X).coefficients(E2X))
sage: fs
x |--> -165/377*cos(4*x) + 240/377*sin(4*x)
sage: fd
x |--> 15/377*(11*cos(x) - 42*sin(x))*e^(-2*x)


The plot is therefore trivial :

sage: plot((fs, fd, fs+fd),(x, 0, 2*pi))
Launched png viewer for Graphics object consisting of 3 graphics primitives HTH,

EDIT : If you are interested in the trigonometric expressions manipulation, here's how to transform your expressions :

sage: a, b, c, d = var("a, b, c, d")


The expected form :

sage: Eq=a*cos(x) + b*sin(x) == c*cos(d + x) ; Eq
a*cos(x) + b*sin(x) == c*cos(d + x)


This relation must be true for any value of $\sin x$ and $cos x$ ; it is therefore implicitly a system of (polynomial equations in $\sin x$ and $cos x$, which we can get explicitly :

sage: Sys=[(Eq.lhs()-Eq.rhs()).trig_expand().coefficient(u(x),1) for u in (sin, cos)] ; Sys
[c*sin(d) + b, -c*cos(d) + a]


Neither Maxima nor Giac or Fricas can directly solve this system, which we will solve manually. Get a first expression of $c$ :

sage: Sc=Sys.solve(c) ; Sc
c == a/cos(d)


which we use to get (a function of) $d$ : sage: St=(((Sys.subs(Sc)==0)-b)/a).trig_reduce() ; St tan(d) == -b/a

We can now get $c^2$ :

sage: Sc2=(Sc.solve(cos(d))^2).subs(cos(d)^2==1/(1+tan(d)^2)).subs(St).solve(c^2) ; Sc2
[c^2 == a^2 + b^2]


Note that we have to manually introduce $\tan d$. We can now get our candidate solutions :

sage: Sol=[[St.solve(d), u] for u in Sc2.solve(c)] ; Sol
[[d == -arctan(b/a), c == -sqrt(a^2 + b^2)],
[d == -arctan(b/a), c == sqrt(a^2 + b^2)]]


Since we have used a solution of $c^2$ to get a solution for $c$, we may have introduced a spurious solution. Indeed, the first solution doesn't check :

sage: bool(Eq.subs(Sol).simplify_full().canonicalize_radical())
False


but the second one does : sage: bool(Eq.subs(Sol).simplify_full().canonicalize_radical()) True

Note that Mathematica

• also needs a manual explicitation of the implicit system,

• expresses its results in an awkward form, not easily simplifiable :

i. e. :

 sage: mathematica.Solve([(Eq.lhs()-Eq.rhs()).trig_expand().coefficient(u(x),1)==0 for u in (sin, cos)], [c, d])
{{c -> -(a^2/Sqrt[a^2 + b^2]) - b^2/Sqrt[a^2 + b^2],
d -> ConditionalExpression[ArcTan[-(a/Sqrt[a^2 + b^2]),
b/Sqrt[a^2 + b^2]] + 2*Pi*C, Element[C, Integers]]},
{c -> a^2/Sqrt[a^2 + b^2] + b^2/Sqrt[a^2 + b^2],
d -> ConditionalExpression[ArcTan[a/Sqrt[a^2 + b^2],
-(b/Sqrt[a^2 + b^2])] + 2*Pi*C, Element[C, Integers]]}}

sage: mathematica.Simplify(mathematica.Solve([(Eq.lhs()-Eq.rhs()).trig_expand().coefficient(u(x),1)==0 for u in (sin, cos)], [c, d]))
{{c -> -Sqrt[a^2 + b^2], d -> ConditionalExpression[
ArcTan[-(a/Sqrt[a^2 + b^2]), b/Sqrt[a^2 + b^2]] + 2*Pi*C,
Element[C, Integers]]}, {c -> Sqrt[a^2 + b^2],
d -> ConditionalExpression[ArcTan[a/Sqrt[a^2 + b^2],
-(b/Sqrt[a^2 + b^2])] + 2*Pi*C, Element[C, Integers]]}}

sage: mathematica.FullSimplify(mathematica.Solve([(Eq.lhs()-Eq.rhs()).trig_expand().coefficient(u(x),1)==0 for u in (sin, cos)], [c, d]))
{{c -> -Sqrt[a^2 + b^2], d -> ConditionalExpression[
2*Pi*C - I*Log[-((a - I*b)/Sqrt[a^2 + b^2])],
Element[C, Integers]]}, {c -> Sqrt[a^2 + b^2],
d -> ConditionalExpression[2*Pi*C - I*Log[(a - I*b)/Sqrt[a^2 + b^2]],
Element[C, Integers]]}}


Sympy can solve the same system, but its solution is ... disputable... :

sage: SS = solve(Sys, [c, d], algorithm="sympy") ; SS
[{c: -(a^2 + b^2 - sqrt(a^2 + b^2)*a)/(a - sqrt(a^2 + b^2)),
d: 2*arctan((a - sqrt(a^2 + b^2))/b)},
{c: -(a^2 + b^2 + sqrt(a^2 + b^2)*a)/(a + sqrt(a^2 + b^2)),
d: 2*arctan((a + sqrt(a^2 + b^2))/b)}]


One can note that :

sage: ((SS[c].numerator()*(a^2 + b^2 + sqrt(a^2 + b^2)*a)).expand()/((SS[c].denominator()*(a^2 + b^2 + sqrt(a^2 + b^2)*a)).expand())).factor()
sqrt(a^2 + b^2)
sage: ((SS[c].numerator()*(sqrt(a^2 + b^2)*a-(a^2 + b^2 ))).expand()/((SS[c].denominator()*(sqrt(a^2 + b^2)*a-(a^2 + b^2 ))).expand())).simplify_full()
-sqrt(a^2 + b^2)


and

sage: tan(SS[d]).trig_expand().expand().canonicalize_radical()
-b/a
-b/a


The solutions given by Sympy are therefore identical to those manually obtained ; the same limitations apply. However, their properties are way less obvious, and their transformation to "easier" forms quite unintuitive...

EDIT 2: Pulling it all together :

f = function('f')
de = diff(f(x),x,2) + 4*diff(f(x),x) + 5*f(x) == 15*cos(4*x)
PM = desolve(de, f(x), ivar=x, ics=(0,0,0))
PM0, PM1 = ((u*(e^(-2*x))^u) for u in PM.coefficients(e^(-2*x)))
var('a, b, c, d')
Eq0 = PM0 - a*cos(4*x + b).trig_expand()
S0 = solve([Eq0.coefficient(u(4*x), 1) for u in (sin, cos)],[a, b],
algorithm = "sympy")
check0 = all([Eq0.subs(s).simplify_full() == 0 for s in S0])
fs = (a*cos(4*x + b)).subs(S0).function(x)
Eq1 = (PM1/e^(-2*x) - c*cos(x+d).trig_expand())
S1=solve([Eq1.coefficient(u(x), 1) for u in (sin, cos)], (c, d), algorithm="sympy")
check1 = all([Eq1.subs(s).simplify_full() == 0 for s in S1])
fd = ((c*cos(x + d)).subs(S1)*e^(-2*x)).function(x)
ft = fs + fd


The result can be seen as :

plot((fs, fd, ft), (0, 2*pi), legend_label=["fs", "fd", "ft"], axes=False, frame=True)


HTH,

more

Thank you for your in depth guide, but my teacher wants the X-spring(t) and X(t) in Ccos(u+d) form, not Acos(u)+Bsin(u) form. I should've specified in my original question, but I was looking for how to convert the original answer I got into the form of the answer that my teacher has

1

I should've specified in my original question,

Indeed.

My edit showed you two ways Sage can be used to compute the relevant values, leaving only the substitution itself. Do you need help with that ? See my second edit.