1

I'm going to solve a ODE system on the form: $$\dot x(t) = f(x(t),u(t))$$

Where an example of the system migth look like:

$$(\dot x_1(t) ,\dot x_2(t) ,\dot x_3(t))= a x_1(t) + b x_2 (t) + c x_3 (t) + d u(t) e x_2 (t) + f u(t) g x_1(t) + h sin(x_3(t))$$

The parameters $a,b,c,d,e,f,g,h$ are known and $x_1,x_2,x_3,\dot x_1,\dot x_2,\dot x_3$ are known for $t = 0$ and $u$ are known $\forall t$. Also $u$ is constant $\forall t$.

So is there a way to to find $x_1,x_2,x_3$ by using linear algebra? I just got started with Armadillo and C++, then I realize that Armadillo does not have and ODE solver. But Armadillo is optimized for linear algebra.

euraad
  • 2,934
  • The example system seems to be written incorrectly. On the left, you should have $(x_1'(t), x_2'(t), x_3'(t))$. You could implement your own ODE solver (using a method like RK4), or alternatively you could probably find a C++ library that provides an ODE solver that you can use. Here's a thread with some library suggestions: https://stackoverflow.com/questions/7622286/is-there-a-c-library-for-ordinary-differential-equation-ode-solvers – littleO Dec 17 '18 at 00:12
  • I'm used to write the x_1 with dots above. Is Euler better? – euraad Dec 17 '18 at 00:14
  • The issue is not with the dots, it's that on the left you've written a sum $x_1'(t) + x_2'(t) + x_3'(t)$ rather than a vector $(x_1'(t), x_2'(t), x_3'(t))$. If I understand correctly, $x$ is a vector-valued function with component functions $x_1, x_2$, and $x_3$. The dot notation is fine. – littleO Dec 17 '18 at 00:16
  • Okej! Maybe I should change. Who is fastest to compute, Euler or RK4? – euraad Dec 17 '18 at 00:19
  • I think that's a question for stackoverflow – Andrei Dec 17 '18 at 00:30
  • 2
    Euler's method is easier to implement, but RK4 is much more accurate. I suspect that RK4 will be faster because you will be able to use a larger stepsize to achieve the same accuracy. There are more sophisticated methods such as RK45 that use adaptive stepsizes, and also implicit methods that are useful for "stiff" ODEs (for which explicit methods like Euler or RK4 might fail). You can implement RK4 pretty easily, but, using an ODE library is probably a good way to go. – littleO Dec 17 '18 at 00:31
  • @littleO Thank you! I will look for a library. – euraad Dec 17 '18 at 00:32
  • @littleO By the way. It need to be fast. I'm going to simulate the same ODE with different $u$ at 256 times to look which $u$ gives most value for the output. I hope C++ will do that under 0.1 seconds. I'm only using small amout of data. – euraad Dec 17 '18 at 00:38
  • I've discussed the issue of accuracy vs. computational effort in the answer to Quick question ... about Runge-Kutta (see also motivation of Runge-Kutta). In short, if you halve the step size in an order $p$ method the error reduces by $2^p$. To get the same error reduction in the Euler method, you have to reduce the step size by a factor $2^p$. On top of that, for any reasonable accuracy larger error range, the higher order method will still take less steps than Euler. – Lutz Lehmann Dec 17 '18 at 00:56
  • Look also for literature on "optimal control". This type of optimization problem for a control function and how to look most efficiently for numerical solutions has some decades of research (and expressive jargon like "bang-bang solutions") behind it. – Lutz Lehmann Dec 17 '18 at 00:59
  • @LutzL Are you sure? What I can see here, somebody claims that Euler is faster that RK45 https://math.stackexchange.com/questions/173718/fastest-numeric-method-for-ode – euraad Dec 17 '18 at 10:38
  • @DanielMårtensson : For a fixed step-size, yes, Euler requires only one function evaluation, higher order methods have $s=p+o(p)$ stages with one function evaluation each. However, for a comparable error level $ε$ Euler requires $O(1/ε)$ steps while the higher order method only requires $O(1/ε^{1/p})$. As I wrote in the links, to get to $ε=10^{-12}$ on the interval $[0,1]$, Euler fails as that would require $10^{12}$ steps where the floating point noise accumulates to the level of about $10^{-4}$. RK4 get there with step size $h=10^{-3}$ with only $4\cdot 10^3=4000$ function evaluations. – Lutz Lehmann Dec 17 '18 at 10:50
  • @DanielMårtensson : For a fixed step-size, yes, Euler requires only one function evaluation, higher order methods have $s=p+o(p)$ stages with one function evaluation each. However, for a comparable error level $ε$ Euler requires $O(1/ε)$ steps while the higher order method only requires $O(1/ε^{1/p})$. As I wrote in the links, to get to $ε=10^{-12}$ on the interval $[0,1]$, Euler fails as that would require $10^{12}$ steps where the floating point noise accumulates to the level of about $10^{-4}$. RK4 get there with step size $h=10^{-3}$ with only $4\cdot 10^3=4000$ function evaluations. – Lutz Lehmann Dec 17 '18 at 10:50
  • However, in the method of lines one has other concerns to consider, as interpolation has the complexity of an Euler step, a long RK4 step with interpolations in between will be slower than pure Euler method. However, if you are not interested in the result of every Euler step, if you have to reduce the step size to get to some required accuracy and thus group Euler steps, then the higher order method becomes more advantageous. -- The space discretization results in a lower bound on accuracy. If that error is $Δx^2$, it makes no sense to drive the error of the time integration much lower. – Lutz Lehmann Dec 17 '18 at 10:58
  • @LutzL Hmmm...I think I could use Euler because I don't care about precision and adaptive step size. The reason why is I only want to simulate an ODE very very fast. As fast as possible. – euraad Dec 17 '18 at 10:59
  • @LutzL I will simulate and second order ODE at 1:10 steps for 256 times. I hope I can do that in C++ for under 0.1 seconds. – euraad Dec 17 '18 at 11:00
  • Well, if you are not afraid of nonsense results, star systems imploding or exploding etc., go ahead. If there are conserved quantities, it helps a little, but not too much, to force them to be constant. This is discussed for instance in the paper/book by Hairer et al. on geometric integration. – Lutz Lehmann Dec 17 '18 at 11:02
  • @LutzL Yes. I know. If where are another way to find the best input $u$ for the ODE, let me know. :) – euraad Dec 17 '18 at 11:06
  • Please edit your question so that the problem statement is complete. Obviously you have some quality function that evaluates on the final state of the integration. Perhaps give also a more complete mock-up problem, for example the inverted pendulum you already used. Then give some keywords on the search method you employ, is it some kind of genetic algorithm, ant hill, tabu search,... or is some kind of gradient search or other derivative based method involved,...? It would also be of interest where you expect the large dimensional linear algebra to occur. – Lutz Lehmann Dec 17 '18 at 11:18

0 Answers0