I have a nonlinear model ${\bf{\dot x}} = f\left( {\bf{x}} \right) $ which is represented as follows:

$$\begin{array}{l}{{\dot x}_1} = {x_2}\\{{\dot x}_2} = {x_3} - \alpha {x_4}{x_2} + \frac{\alpha }{T}{x_1}\\{{\dot x}_3} = \left( { - \alpha - \frac{\alpha }{T}} \right){x_2}\\{\dot x_4} = - \frac{1}{T}{x_4}\end{array} % MathType!MTEF!2!1!+- % feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn % hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr % 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9 % vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x % fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGceaqabeaaceWG4b % GbaiaadaWgaaWcbaGaaGymaaqabaGccqGH9aqpcaWG4bWaaSbaaSqa % aiaaikdaaeqaaaGcbaGabmiEayaacaWaaSbaaSqaaiaaikdaaeqaaO % Gaeyypa0JaamiEamaaBaaaleaacaaIZaaabeaakiabgkHiTiabeg7a % HjaadIhadaWgaaWcbaGaaGinaaqabaGccaWG4bWaaSbaaSqaaiaaik % daaeqaaOGaey4kaSYaaSaaaeaacqaHXoqyaeaacaWGubaaaiaadIha % daWgaaWcbaGaaGymaaqabaaakeaaceWG4bGbaiaadaWgaaWcbaGaaG % 4maaqabaGccqGH9aqpdaqadaqaaiabgkHiTiabeg7aHjabgkHiTmaa % laaabaGaeqySdegabaGaamivaaaaaiaawIcacaGLPaaacaWG4bWaaS % baaSqaaiaaikdaaeqaaaGcbaGaamiEamaaBaaaleaacaaI0aaabeaa % kiabg2da9iabgkHiTmaalaaabaGaaGymaaqaaiaadsfaaaGaamiEam % aaBaaaleaacaaI0aaabeaaaaaa!5F81! $$

where $\alpha, T$ are constants.
I want to calculate ${{\bf{x}}^{\left( n \right)}} = f\left( {\bf{x}} \right)$.
I am new to Mathematica, so any reference or help is appreciated. Thanks.

  • $\begingroup$ Please explain how discrete steps are obtained from these differential equations. Also, should the left side of the fourth equation be a first derivative? $\endgroup$
    – bbgodfrey
    Commented Mar 7, 2022 at 6:02
  • $\begingroup$ Hello, I modified the question, sorry for the unclarity, I am interested in calculating the nth derivative of the system equations. $\endgroup$
    – TanjiroLL
    Commented Mar 7, 2022 at 7:43
  • $\begingroup$ This model is linear; you can see this by solving equation 4 as $x_4(t)=Ae^{-t/T}$ and substituting back into equation 2. $\endgroup$
    – Jojo
    Commented Mar 7, 2022 at 16:09

3 Answers 3


I would like to calculate the system state after n steps

To do that, use OutputResponse. This gives the output of the system. You decide on how long to look at the output. Also your system is continuous not discrete.

But to do this, you need to have $u$ added to your equations (the external forcing function). Otherwise the system will not move. (since no initial conditions are given, system assumed at rest at time $t=0$).

Right now you do not have any. Then you can tell Mathematica to solve the system and show the output with some specific input to use for u. Like this (notice, added $u$ for first equation. You can change that)


Mathematica graphics

Now you can move it. Used UnitStep for $u$ here. You can change that also.

Plot[output, {t, 0, 5}, AxesLabel -> {"time", "state x"}, 
 BaseStyle -> 12, GridLines -> Automatic, 
 GridLinesStyle -> LightGray]

Mathematica graphics

Each curve above is the output of each state variable.


I am interested in calculating the derivatives, not getting system response

If you want to just take derivatives, then you do not need state control. Just use D operator

ClearAll[x1, x2, x3, x4, t,n]
sys = {x2[t], x3[t] - alpha*x4[t]*x2[t] + alpha/T*x1[t], (-alpha - alpha/T)*x2[t], -1/T*x4[t]}

Now take as many derivatives as you want

D[sys, t]

Mathematica graphics

D[sys, {t, 2}]

Mathematica graphics

D[sys, {t, n}]

Mathematica graphics

  • $\begingroup$ Thank you for your help. First point: Yes, the given equations are continuous, because I want to calculate the nth derivative, then I will do the discretization. Second point: for input, u is constant, that's why I didn't added in the equations, Third points: I am interested in calculating the derivatives, not getting system response, because I am controlling the system on matlab using NMPC controller $\endgroup$
    – TanjiroLL
    Commented Mar 7, 2022 at 7:36
  • $\begingroup$ @coder00 fyi, added just taking derivatives. $\endgroup$
    – Nasser
    Commented Mar 7, 2022 at 7:54

Suppose that we have $x'=f(x), x^{(n-1)}=F_{n-1}(x)$, then $x^{(n)}=x'.\nabla F_{n-1}=f.\nabla F_{n-1}$, where $\nabla_i=\partial_{x_i}$. therefore, with given x'=f[x] we can compute any derivative $x^{(n)}=F_n(x)$ as follows:

f = {x[2], 
   x[3] - a x[4] x[2] + a/T x[1], -a (1 + 1/T) x[2], -1/T x[4]};
q[y_] := Sum[f[[i]] D[y, x[i]], {i, 4}];

Second derivative

d2x=q[f] // Simplify

(*Out[]= {(a x[1])/T + x[3] - a x[2] x[4], 
 a (-x[2] + ((-a x[1] + x[2]) x[4])/T - x[3] x[4] + a x[2] x[4]^2), (
 a (1 + T) (-a x[1] - T x[3] + a T x[2] x[4]))/T^2, x[4]/T^2}*)

Third derivative

 d3x = q[q[f]] // Simplify

(*Out[]= {a (-x[2] + ((-a x[1] + x[2]) x[4])/T - x[3] x[4] + 
    a x[2] x[4]^2), -((
  a (T^2 x[3] + x[2] x[4] - 2 T x[3] x[4] + 
     a^2 T x[4]^2 (-x[1] + T x[2] x[4]) + 
     a (-2 x[1] x[4] - T^2 x[4] (2 x[2] + x[3] x[4]) + 
        T (x[1] + 3 x[2] x[4]^2))))/T^2), -((
  a^2 (1 + T) ((-a x[1] + x[2]) x[4] - 
     T (x[2] + x[3] x[4] - a x[2] x[4]^2)))/T^2), -(x[4]/T^3)}*)

Fourth derivative

d4x = q[q[q[f]]] // Simplify

(*Out[]= {-((
  a (T^2 x[3] + x[2] x[4] - 2 T x[3] x[4] + 
     a^2 T x[4]^2 (-x[1] + T x[2] x[4]) + 
     a (-2 x[1] x[4] - T^2 x[4] (2 x[2] + x[3] x[4]) + 
        T (x[1] + 3 x[2] x[4]^2))))/T^2), (1/(T^3))
 a ((x[2] - 3 T x[3]) x[4] + a^3 T^2 x[4]^3 (-x[1] + T x[2] x[4]) + 
    a (-3 x[1] x[4] + 7 T x[2] x[4]^2 + T^3 (x[2] + 2 x[3] x[4]) - 
       T^2 x[4] (4 x[2] + 5 x[3] x[4])) - 
    a^2 T x[4] (5 x[1] x[4] + T^2 x[4] (3 x[2] + x[3] x[4]) - 
       2 T (x[1] + 3 x[2] x[4]^2))), (1/(T^3))
 a^2 (1 + T) (T^2 x[3] + x[2] x[4] - 2 T x[3] x[4] + 
    a^2 T x[4]^2 (-x[1] + T x[2] x[4]) + 
    a (-2 x[1] x[4] - T^2 x[4] (2 x[2] + x[3] x[4]) + 
       T (x[1] + 3 x[2] x[4]^2))), x[4]/T^4}*)

Technically we can organize this series with Table,Do or NestList.


The r.h.s and the equations expressed as rules.

f = {x2[t], (a x1[t])/T + x3[t] - α x2[t] x4[t], (-(a/T) - α) x2[t], -(x4[t]/T)}
rules = Thread[D[{x1[t], x2[t], x3[t], x4[t]}, t] -> f]

Take the first derivative of $f$ and substitute the rules to get ${{\bf{x}}^{\left( 2 \right)}}$.

D[f, t] /. rules

To do the same thing for the third derivative, repeat it for the above result.

D[%, t] /. rules

This is a Nest operation. To get to any derivative, say 3, do

Nest[D[#, t] /. rules &, f, 3-1]

$$\left\{\frac{\alpha \left(\text{x2}(t) \left(\alpha T \text{x4}(t)^2+\text{x4}(t)-T\right)-\text{x4}(t) (a \text{x1}(t)+T \text{x3}(t))\right)}{T},\frac{\alpha \left(a \text{x1}(t) (\text{x4}(t) (\alpha T \text{x4}(t)+2)-T)-\text{x2}(t) \text{x4}(t) \left(\alpha T \text{x4}(t) (\alpha T \text{x4}(t)+3)-2 \alpha T^2+1\right)+T \text{x3}(t) (\text{x4}(t) (\alpha T \text{x4}(t)+2)-T)\right)}{T^2},\frac{\alpha (a+\alpha T) \left(\text{x4}(t) (a \text{x1}(t)+T \text{x3}(t))-\text{x2}(t) \left(\alpha T \text{x4}(t)^2+\text{x4}(t)-T\right)\right)}{T^2},-\frac{\text{x4}(t)}{T^3}\right\}$$

If you need the intermediate results use NestList.

These will give any derivative as a function of the original states ${{\bf{x}}^{\left( n \right)}}=g\left(\bf{x}\right)$. If the system has inputs, the derivatives of these inputs will appear in $g$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.