Linearization of differential equation

Questions
1

Hello,

I would like to linearize a differential equation around a equilibrium position. The description of the steps that I have carried out are :

Equation

eq = 1/2 g l m Cos[?[t]] == J (?^??)[t]

Equation around the equilibrium position can be re-written as :

eqAe = Replace[eq, ?[ t] -> (?e &) + ??[t], ?]

Change of variables

eqLi = Replace[eq, {D[??[t], {t, 2}] -> ??pp, D[??[t], t] -> ??p, ??[t] -> ??}, ?]

Linearization

  eqLi2 = Series[eqLi[[1]], {{??pp, ??p,  ??}, {0, 0, 0}, 1}] == 
          Series[eqLi[[2]], {{??pp, ??p, ??}, {0, 0, 0}, 1}]

Change of variables

eqLi3 = Replace[eq, ??pp -> {D[??[t], {t, 2}], ??p -> D[??[t], t], ?? -> ??[t]}, ?]

Currently, I try to transform a code that I have done on Maple and which works to Mathematica. However, I have still some blocking points.

Can anyone help me to obtain to correct my code ?
I send the .nb in copy

Thanks a lot for your help

2

I would try something like this:

eqAe = eq /. \[Theta] -> (\[Theta]e + \[CapitalDelta]\[Theta][#] &)
repl = {\[CapitalDelta]\[Theta]''[
     t] -> \[CapitalDelta]\[Theta]pp, \[CapitalDelta]\[Theta]'[
     t] -> \[CapitalDelta]\[Theta]p, \[CapitalDelta]\[Theta][
     t] -> \[CapitalDelta]\[Theta]};
EqLi = eqAe /. repl
Normal@Series[
   EqLi, {\[CapitalDelta]\[Theta]pp, 0, 1}, {\[CapitalDelta]\[Theta]p,
     0, 1}, {\[CapitalDelta]\[Theta], 0, 1}] /. Map[Reverse, repl]
3

Usually you work with Jacobian matrix to find the equilibrium point of system via linearization. You will have much simpler set of equations to deal with. y1 = \theta and y2=\theta'. [y1',y2']^{T} = f([y1,y2]^{T}). The Jacobian is hence:

D[{y2, (g l m Cos[y1])/(2 J)}, {{y1, y2}}]
====> {{0, 1}, {-((g l m Sin[y1])/(2 J)), 0}}

you find the theta at equlibrium by

Det[%] ====> (g l m Sin[y1])/(2 J)

Therefore the equilibrium occurs at \sin{y1} ==0 or \sin{\theta} == 0

You may refer to 2.Stability of this article and Vector Analysis in Mathematica

4

Thank you Gianluca Gorni.
You perfectly answer to my question.
However, I have still a slight question.
In fact, I would have liked to linearize the equations around the equilibrium position minus the equation at the equilibrium position. I have attached a .nb in order to give you more precisions.
The point is : Is it possible to make soustraction between 2 equations ? How can we code equation 3 = equation 2 - equation 1?
May you have a look on my file to see how it is possible to correct the soustraction between 2 equations.
Thank you for your help.

5

You can combine equations using rules:

equationCombinationRules = {Equal[a_, b_] + Equal[c_, d_] :> 
    Equal[a + c, b + d],
   a_?NumericQ*Equal[c_, d_] :> Equal[a*c, a*d]};

-(1/2 g l m Cos[\[Theta]e] == 0) + (1/2 g l m Cos[\[Theta]e] - 
     1/2 g l m \[CapitalDelta]\[Theta] Sin[\[Theta]e] == 
    J \[CapitalDelta]\[Theta]pp) //. equationCombinationRules
6

OK. Thanks a lot for your help.
But is not there something simpler ?
I find this code interesting but a bit coimplex just to make a soustraction between two equations

I have almost a solution with this code :
e1 = 3 a + e == b;
e2 = 6 c == d;
Thread[e1 + e2, Equal]

However with minus, it doesn’t work … :slight_smile:

7

Yes, it is a bit complicated. If you need it a lot, you can load the following definitions in your initialization code:

Unprotect[Equal];
Equal /: Equal[a_, b_] + Equal[c_, d_] := Equal[a + c, b + d];
Equal /: a_?NumericQ*Equal[c_, d_] := Equal[a*c, a*d];
Protect[Equal];

After that, you can simply type

(a == b) - 2 (c == d)

and get a - 2 c == b - 2 d as desired.

I think Wolfram should provide this behaviour as default, if there is no serious reason not to do it.

8

Hello,

Thank you for your help.

First, concerning the substraction of two equations, I found on the net on more simple way (to my point of view) :

 `Inner[Subtract, eq\[Psi]Ae, eq\[Psi]E, Equal];`

Secondly, I have still a little concerning the linearization that I have noticed when I have employed my method on a more complex example (see in the file attached).

I find this for the equation of psi

\Delta \psi '(t) \left(\Delta \theta (t) \Delta \theta '(t) \left(-2 A \sin (\text{$\theta $e}) \cos '(\text{$\theta $e})-2 A \cos (\text{$\theta $e}) \sin '(\text{$\theta $e})+2 C \sin (\text{$\theta $e}) \cos '(\text{$\theta $e})+2 C \cos (\text{$\theta $e}) \sin '(\text{$\theta $e})\right)+\Delta \theta '(t) (2 C \sin (\text{$\theta $e}) \cos (\text{$\theta $e})-2 A \sin (\text{$\theta $e}) \cos (\text{$\theta $e}))\right)+\Delta \psi ''(t) \left(\Delta \theta (t) \left(2 A \cos (\text{$\theta $e}) \cos '(\text{$\theta $e})-2 C \cos (\text{$\theta $e}) \cos '(\text{$\theta $e})\right)+A \cos (\text{$\theta $e})^2-C \cos (\text{$\theta $e})^2+C\right)+C \Omega \Delta \theta (t) \cos '(\text{$\theta $e}) \Delta \theta '(t)+C \Omega \cos (\text{$\theta $e}) \Delta \theta '(t)=\Delta \psi (t) \left(g L m_3 \cos (\text{$\theta $e}) \sin '(\text{$\psi $e})+g L m_3 \Delta \theta (t) \cos '(\text{$\theta $e}) \sin '(\text{$\psi $e})\right)+g L m_3 \Delta \theta (t) \sin (\text{$\psi $e}) \cos '(\text{$\theta $e})

and for the equation theta :

A \Delta \theta ''(t)+\Delta \psi '(t) \left(-C \Omega \cos (\text{$\theta $e})-C \Omega \Delta \theta (t) \cos '(\text{$\theta $e})\right)=\Delta \psi (t) \left(g L m_3 \sin (\text{$\theta $e}) \cos '(\text{$\psi $e})+g L m_3 \Delta \theta (t) \sin '(\text{$\theta $e}) \cos '(\text{$\psi $e})\right)+g L m_3 \Delta \theta (t) \cos (\text{$\psi $e}) \sin '(\text{$\theta $e})

Instead of

psi equation :

C \left( {\frac {\rm d}{{\rm d}t}}\epsilon2 \left( t \right) \right) \Omega\,\cos \left( {\it theta\_e} \right) + \left( A \left( \cos \left( {\it theta\_e} \right) \right) ^{2}-C \left( \cos \left( { \it theta\_e} \right) \right) ^{2}+C \right) {\frac {{\rm d}^{2}}{ {\rm d}{t}^{2}}}\epsilon1 \left( t \right) =\epsilon1 \left( t \right) Lgm_{{3}}\cos \left( {\it psi\_e} \right) \cos \left( {\it theta\_e} \right) -\epsilon2 \left( t \right) Lgm_{{3}}\sin \left( { \it psi\_e} \right) \sin \left( {\it theta\_e} \right)

for theta equation

-C \left( {\frac {\rm d}{{\rm d}t}}\epsilon1 \left( t \right) \right) \Omega\,\cos \left( {\it theta\_e} \right) +A{\frac {{\rm d}^ {2}}{{\rm d}{t}^{2}}}\epsilon2 \left( t \right) =-\epsilon1 \left( t \right) Lgm_{{3}}\sin \left( {\it psi\_e} \right) \sin \left( {\it theta\_e} \right) +\epsilon2 \left( t \right) Lgm_{{3}}\cos \left( { \it psi\_e} \right) \cos \left( {\it theta\_e} \right)

So, it seems that i have some terms with second orders which are still in my equations after linearization at the first order.

Do you have ideas how I can correct my linearization step ?

Thanks a lot of your help

P.S: Sorry if my post is not easy to read, I have still some difficulties to enter Mathematica inputs in my posts. I’m interesting also how I can enter Mathematica inputs in a post.

9

Sorry, I have no expertise in linearizing differential equations.