ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 03 Jun 2020 23:27:14 +0200Solving 2nd order ode with desolve_system_rk4https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/ Dear All,
I have the following nonlinear ode to solve numerically:
$$\ddot q(t) = -\dot q(t)^2 - q(t)$$.
Since desolve_system_rk4 cannot solve 2nd order ode I make this ode a system of two 1st order odes in the following way:
$$q'(t) = qdot(t)$$
$$qdot'(t) = -qdot(t)^2 - q(t)$$
However desolve_system_rk4 cannot still solve this 2 odes numerically. It raises an error. Any help is appreciated.
Here is the code I wrote:
t = var('t')
q = function('q')(t)
qdot = function('qdot')(t)
# Equations of motion
eom1 = q.diff(t) == qdot(t)
eom2 = qdot.diff(t) == -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], [q, qdot], ivar = t , ics = [q0, qdot0])Mon, 01 Jun 2020 10:58:23 +0200https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/Answer by dsejas for <p>Dear All,</p>
<p>I have the following nonlinear ode to solve numerically:</p>
<p>$$\ddot q(t) = -\dot q(t)^2 - q(t)$$.</p>
<p>Since desolve_system_rk4 cannot solve 2nd order ode I make this ode a system of two 1st order odes in the following way:</p>
<p>$$q'(t) = qdot(t)$$
$$qdot'(t) = -qdot(t)^2 - q(t)$$</p>
<p>However desolve_system_rk4 cannot still solve this 2 odes numerically. It raises an error. Any help is appreciated.</p>
<p>Here is the code I wrote:</p>
<pre><code>t = var('t')
q = function('q')(t)
qdot = function('qdot')(t)
# Equations of motion
eom1 = q.diff(t) == qdot(t)
eom2 = qdot.diff(t) == -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], [q, qdot], ivar = t , ics = [q0, qdot0])
</code></pre>
https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?answer=51692#post-id-51692Hello, @F. Semih Dündar! I can see a couple of problems with your code. First of all, there is no need to define $q$ and $qdot$ as formal functions, since that is implicitly assumed. Second, you only need to pass the RHS (right hand side) of the differential equations. Finally, the `ics` parameter must be of the form `[t0, q0, qdot0]`.
With those considerations, your code should look like this:
var('t, qdot, q') # independent variable and functions
# Equations of motion
eom1 = qdot
eom2 = -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], ivar=t, vars=[q, qdot], ics = [0, q0, qdot0], end_points=1.8, step=0.01)
S = [(t, q) for t, q, qdot in sol]
list_plot(S, plotjoined=True)
I think a few additional clarifications are in order. As mentioned, you only need the RHS of the DEs; that is why I removed part of $oem1$ and $oem2$ (specifically, the LHS and the `==`). The `vars` parameter implicitly indicate that your system is of the form
$$q' = oem1$$
$$qdot' = oem2$$
so this code defines the same system you intended.
I have used the `end_points=1.8` parameter in order to solve in the interval $[0, 1.8]$---after $t\approx1.8$, the value of $q$ is negative, which I assumed not to be desirable, although I don't know the context in which you are trying to solve this (so maybe you would like to modify that parameter if needed).
On the other hand, when I plotted the solution, due to the small domain, I obtained a pointy curve, so I used the `step=0.01` parameter, in order to decrease the step length, and obtain a more refined solution.
The last two lines are intended to obtain a plot of $q$ alone, without including $qdot$. (I assume they are clear to you, but feel free to ask otherwise.)
Here is the resulting image:
![image description](/upfiles/1591059860228864.png)
I hope this helps!Tue, 02 Jun 2020 03:04:39 +0200https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?answer=51692#post-id-51692Comment by aszanti for <p>Hello, <a href="/users/27954/f-semih-dundar/">@F. Semih Dündar</a>! I can see a couple of problems with your code. First of all, there is no need to define $q$ and $qdot$ as formal functions, since that is implicitly assumed. Second, you only need to pass the RHS (right hand side) of the differential equations. Finally, the <code>ics</code> parameter must be of the form <code>[t0, q0, qdot0]</code>.</p>
<p>With those considerations, your code should look like this:</p>
<pre><code>var('t, qdot, q') # independent variable and functions
# Equations of motion
eom1 = qdot
eom2 = -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], ivar=t, vars=[q, qdot], ics = [0, q0, qdot0], end_points=1.8, step=0.01)
S = [(t, q) for t, q, qdot in sol]
list_plot(S, plotjoined=True)
</code></pre>
<p>I think a few additional clarifications are in order. As mentioned, you only need the RHS of the DEs; that is why I removed part of $oem1$ and $oem2$ (specifically, the LHS and the <code>==</code>). The <code>vars</code> parameter implicitly indicate that your system is of the form</p>
<p>$$q' = oem1$$
$$qdot' = oem2$$</p>
<p>so this code defines the same system you intended.</p>
<p>I have used the <code>end_points=1.8</code> parameter in order to solve in the interval $[0, 1.8]$---after $t\approx1.8$, the value of $q$ is negative, which I assumed not to be desirable, although I don't know the context in which you are trying to solve this (so maybe you would like to modify that parameter if needed).</p>
<p>On the other hand, when I plotted the solution, due to the small domain, I obtained a pointy curve, so I used the <code>step=0.01</code> parameter, in order to decrease the step length, and obtain a more refined solution.</p>
<p>The last two lines are intended to obtain a plot of $q$ alone, without including $qdot$. (I assume they are clear to you, but feel free to ask otherwise.)</p>
<p>Here is the resulting image:</p>
<p><img alt="image description" src="/upfiles/1591059860228864.png"></p>
<p>I hope this helps!</p>
https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51744#post-id-51744Hello dsejas! I understand now, thank you very much for this explanation.Wed, 03 Jun 2020 23:27:14 +0200https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51744#post-id-51744Comment by F. Semih Dündar for <p>Hello, <a href="/users/27954/f-semih-dundar/">@F. Semih Dündar</a>! I can see a couple of problems with your code. First of all, there is no need to define $q$ and $qdot$ as formal functions, since that is implicitly assumed. Second, you only need to pass the RHS (right hand side) of the differential equations. Finally, the <code>ics</code> parameter must be of the form <code>[t0, q0, qdot0]</code>.</p>
<p>With those considerations, your code should look like this:</p>
<pre><code>var('t, qdot, q') # independent variable and functions
# Equations of motion
eom1 = qdot
eom2 = -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], ivar=t, vars=[q, qdot], ics = [0, q0, qdot0], end_points=1.8, step=0.01)
S = [(t, q) for t, q, qdot in sol]
list_plot(S, plotjoined=True)
</code></pre>
<p>I think a few additional clarifications are in order. As mentioned, you only need the RHS of the DEs; that is why I removed part of $oem1$ and $oem2$ (specifically, the LHS and the <code>==</code>). The <code>vars</code> parameter implicitly indicate that your system is of the form</p>
<p>$$q' = oem1$$
$$qdot' = oem2$$</p>
<p>so this code defines the same system you intended.</p>
<p>I have used the <code>end_points=1.8</code> parameter in order to solve in the interval $[0, 1.8]$---after $t\approx1.8$, the value of $q$ is negative, which I assumed not to be desirable, although I don't know the context in which you are trying to solve this (so maybe you would like to modify that parameter if needed).</p>
<p>On the other hand, when I plotted the solution, due to the small domain, I obtained a pointy curve, so I used the <code>step=0.01</code> parameter, in order to decrease the step length, and obtain a more refined solution.</p>
<p>The last two lines are intended to obtain a plot of $q$ alone, without including $qdot$. (I assume they are clear to you, but feel free to ask otherwise.)</p>
<p>Here is the resulting image:</p>
<p><img alt="image description" src="/upfiles/1591059860228864.png"></p>
<p>I hope this helps!</p>
https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51694#post-id-51694Thanks for your help and time. It really helped me.Tue, 02 Jun 2020 07:51:44 +0200https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51694#post-id-51694Comment by aszanti for <p>Hello, <a href="/users/27954/f-semih-dundar/">@F. Semih Dündar</a>! I can see a couple of problems with your code. First of all, there is no need to define $q$ and $qdot$ as formal functions, since that is implicitly assumed. Second, you only need to pass the RHS (right hand side) of the differential equations. Finally, the <code>ics</code> parameter must be of the form <code>[t0, q0, qdot0]</code>.</p>
<p>With those considerations, your code should look like this:</p>
<pre><code>var('t, qdot, q') # independent variable and functions
# Equations of motion
eom1 = qdot
eom2 = -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], ivar=t, vars=[q, qdot], ics = [0, q0, qdot0], end_points=1.8, step=0.01)
S = [(t, q) for t, q, qdot in sol]
list_plot(S, plotjoined=True)
</code></pre>
<p>I think a few additional clarifications are in order. As mentioned, you only need the RHS of the DEs; that is why I removed part of $oem1$ and $oem2$ (specifically, the LHS and the <code>==</code>). The <code>vars</code> parameter implicitly indicate that your system is of the form</p>
<p>$$q' = oem1$$
$$qdot' = oem2$$</p>
<p>so this code defines the same system you intended.</p>
<p>I have used the <code>end_points=1.8</code> parameter in order to solve in the interval $[0, 1.8]$---after $t\approx1.8$, the value of $q$ is negative, which I assumed not to be desirable, although I don't know the context in which you are trying to solve this (so maybe you would like to modify that parameter if needed).</p>
<p>On the other hand, when I plotted the solution, due to the small domain, I obtained a pointy curve, so I used the <code>step=0.01</code> parameter, in order to decrease the step length, and obtain a more refined solution.</p>
<p>The last two lines are intended to obtain a plot of $q$ alone, without including $qdot$. (I assume they are clear to you, but feel free to ask otherwise.)</p>
<p>Here is the resulting image:</p>
<p><img alt="image description" src="/upfiles/1591059860228864.png"></p>
<p>I hope this helps!</p>
https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51713#post-id-51713Why is "qdot" in the equation eom1? Maybe it should be like F.SemihDündar wrote, i.e. eom1 = qdot (t). In this case, we will receive another solution.Wed, 03 Jun 2020 00:28:11 +0200https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51713#post-id-51713Comment by dsejas for <p>Hello, <a href="/users/27954/f-semih-dundar/">@F. Semih Dündar</a>! I can see a couple of problems with your code. First of all, there is no need to define $q$ and $qdot$ as formal functions, since that is implicitly assumed. Second, you only need to pass the RHS (right hand side) of the differential equations. Finally, the <code>ics</code> parameter must be of the form <code>[t0, q0, qdot0]</code>.</p>
<p>With those considerations, your code should look like this:</p>
<pre><code>var('t, qdot, q') # independent variable and functions
# Equations of motion
eom1 = qdot
eom2 = -(qdot(t))^2 - q(t)
# Initial conditions
q0 = 0
qdot0 = 1
sol = desolve_system_rk4([eom1, eom2], ivar=t, vars=[q, qdot], ics = [0, q0, qdot0], end_points=1.8, step=0.01)
S = [(t, q) for t, q, qdot in sol]
list_plot(S, plotjoined=True)
</code></pre>
<p>I think a few additional clarifications are in order. As mentioned, you only need the RHS of the DEs; that is why I removed part of $oem1$ and $oem2$ (specifically, the LHS and the <code>==</code>). The <code>vars</code> parameter implicitly indicate that your system is of the form</p>
<p>$$q' = oem1$$
$$qdot' = oem2$$</p>
<p>so this code defines the same system you intended.</p>
<p>I have used the <code>end_points=1.8</code> parameter in order to solve in the interval $[0, 1.8]$---after $t\approx1.8$, the value of $q$ is negative, which I assumed not to be desirable, although I don't know the context in which you are trying to solve this (so maybe you would like to modify that parameter if needed).</p>
<p>On the other hand, when I plotted the solution, due to the small domain, I obtained a pointy curve, so I used the <code>step=0.01</code> parameter, in order to decrease the step length, and obtain a more refined solution.</p>
<p>The last two lines are intended to obtain a plot of $q$ alone, without including $qdot$. (I assume they are clear to you, but feel free to ask otherwise.)</p>
<p>Here is the resulting image:</p>
<p><img alt="image description" src="/upfiles/1591059860228864.png"></p>
<p>I hope this helps!</p>
https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51735#post-id-51735Hello, @aszanti! That is not correct. Observe the order in which the first and third arguments for `desolve_system_rk4` are given: the first one is `[oem1, oem2]`, and the second (`vars`) is `[q, qdot]`. As I mentioned in my answer, the syntax of the command then implies that
$$q'=oem1=qdot$$
$$qdot′(t)=oem2=−qdot(t)2−q(t)$$
which is exactly the system that @F. Semih Dündar wants to solve. In particular, `sol` will have the form
$$sol = [t, q(t), qdot(t)],$$
for $t\in[0,1.8]$.
Remember, $oem1$ is not an equation, but the the RHS of an equation. The LHS is implicitly equal to the derivative of the first argument of the `vars` parameter.Wed, 03 Jun 2020 20:05:38 +0200https://ask.sagemath.org/question/51668/solving-2nd-order-ode-with-desolve_system_rk4/?comment=51735#post-id-51735