Ask Your Question

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

asked 2021-04-13 02:05:34 +0100

Presto.creator gravatar image

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


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)

But the answer is actually

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 flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2021-04-13 13:43:30 +0100

Emmanuel Charpentier gravatar image

updated 2021-04-14 07:32:26 +0100

[ 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[0]*E2X^u[1]).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

image description


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[1].solve(c)[0] ; Sc
c == a/cos(d)

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

We can now get $c^2$ :

sage: Sc2=(Sc.solve(cos(d))[0]^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)[0], u] for u in Sc2[0].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[0]).simplify_full().canonicalize_radical())

but the second one does : sage: bool(Eq.subs(Sol[1]).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[1], Element[C[1], 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[1], Element[C[1], 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[1], 
    Element[C[1], 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[1], Element[C[1], 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[1] - I*Log[-((a - I*b)/Sqrt[a^2 + b^2])], 
    Element[C[1], Integers]]}, {c -> Sqrt[a^2 + b^2], 
  d -> ConditionalExpression[2*Pi*C[1] - I*Log[(a - I*b)/Sqrt[a^2 + b^2]], 
    Element[C[1], 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[0][c].numerator()*(a^2 + b^2 + sqrt(a^2 + b^2)*a)).expand()/((SS[0][c].denominator()*(a^2 + b^2 + sqrt(a^2 + b^2)*a)).expand())).factor()
sqrt(a^2 + b^2)
sage: ((SS[1][c].numerator()*(sqrt(a^2 + b^2)*a-(a^2 + b^2 ))).expand()/((SS[1][c].denominator()*(sqrt(a^2 + b^2)*a-(a^2 + b^2 ))).expand())).simplify_full()
-sqrt(a^2 + b^2)


sage: tan(SS[0][d]).trig_expand().expand().canonicalize_radical()
sage: tan(SS[1][d]).trig_expand().expand().canonicalize_radical()

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[0]*(e^(-2*x))^u[1]) 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[0]).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[0])*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)


edit flag offensive delete link 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

Presto.creator gravatar imagePresto.creator ( 2021-04-13 22:53:19 +0100 )edit

I should've specified in my original question,


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.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2021-04-14 07:25:43 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2021-04-13 02:05:34 +0100

Seen: 233 times

Last updated: Apr 14 '21