# Solving differential equations for simulations

Recently I've started a few programming projects which involve simulating physical systems. However, more than most involve a differential equation which needs to be solved. Some examples are:

$$\frac{d^2\theta}{dt^2} = -\frac{g}{l}\sin \theta$$ $$\left(\frac{da}{dt}\times\frac{1}{a}\right)^2 = \frac{8\pi G\rho}{3} - \frac{kc^2}{a^2}$$

The technique I have been doing so far is using Wolfram Alpha to integrate the equation and then plot the result. For some of the equations I have found, WA has struggled, leaving me clueless.

I don't know a huge amount about differential equations beyond simple equations such as $\frac{dy}{dx}=2xy$, so I am unable to integrate the equations by hand.

That was why I was wondering if there is a way to find the solution on the fly, so to speak. For example, I would like to generate a graph of $\theta(t)$ in the first equation without first integrating it.

I mostly use JavaScript and Python, but am willing to use another, free, language which has the tools which will help me complete my projects.

• Have you looked to a free CAS like SAGE or a Matlab clone Octave or many others and using the built-in numerical DEQ solvers?
– Moo
Aug 21, 2017 at 18:01
• @Moo Not yet. I didn't know Octave had anything like that Aug 21, 2017 at 18:03
• It has numerical solvers like Matlab and is free.
– Moo
Aug 21, 2017 at 18:04
• @Moo Nice, I'll have a look Aug 21, 2017 at 18:05
• solving your first equation i got something with Jacobi Amplitude Aug 21, 2017 at 18:05

ODE is your friend, available in:

In MATLAB/Octave, Others solvers ode23, ode15 are intended for specific cases of ODEs, when the default choice of ode45 is not working well, or when the problems include nonlinearities, discontinuities, etc.

There is no difficulty writing a Runge-Kutta solver (RK4), for equations $y'=f(y,t)$ where $f$ is arbitrary. https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

The scalar version is suited to first order ODEs.

For higher order, you turn the equation in a vector form:

$$y'''=f(y'',y',y,t)$$ becomes

$$(y_2,y_1,y_0)'=(f(y_2,y_1,y_0,t),y_2,y_1)$$ which is of the type

$$(\vec y)'=\vec g(\vec y,t).$$

• So I presume $y'' = f(y', y, t)$? Aug 21, 2017 at 21:24
• @BetaDecay: what do you mean ?
– user65203
Aug 21, 2017 at 21:26
• You only show how to use the method on a first order differential equation and a third order. How about a second order? Aug 21, 2017 at 22:43
• @BetaDecay: come on, the generalization to any order is obvious.
– user65203
Aug 21, 2017 at 23:18

With the first one $$\frac{d^{2}\theta(t)}{dt^{2}}=-g\sin{\theta(t)}$$ Multiply by $\frac{d\theta(t)}{dt}$ $$\frac{d\theta(t)}{dt}\frac{d^{2}\theta(t)}{dt^{2}}=-g\frac{d\theta(t)}{dt}\sin{\theta(t)}$$ Using the chain rule $$\frac{d}{dt}\frac{1}{2}\Big(\frac{d\theta(t)}{dt}\Big)^{2}=g\frac{d}{dt}\cos{\theta(t)}$$ Thus $$\Big(\frac{d\theta(t)}{dt}\Big)^{2}=2g\cos{\theta(t)}+c_{1}$$ Hence $$t+c_{2}=\frac{1}{\sqrt{2g}}\int^{\theta}\frac{d\theta'}{\pm\sqrt{\frac{c_{1}}+\cos{\theta'}}}=\frac{2\sqrt{\frac{c_{1}+2g\cos{\theta}}{c_{1}+2g}}F\Big(\frac{\theta}{2}\Big{|}\frac{4g}{c_{1}+1}\Big)}{\pm\sqrt{c_{1}+2g\cos{\theta}}}$$ Where $F(z|k^{2})$ is the incomplete elliptic integral of the first kind. With the second one $$\frac{da}{dt}=\pm\sqrt{\frac{8\pi{G}\rho}{3}a-\frac{kc^{2}}{a}}$$ $$t+c=\int^{a}\frac{da'}{\pm\sqrt{\frac{8\pi{G}\rho}{3}a'-\frac{kc^{2}}{a'}}}=\int^{a}\frac{\sqrt{a'}da'}{\pm\sqrt{\frac{8\pi{G}\rho}{3}(a')^{2}-kc^{2}}}$$ Let $a'=\sqrt{\frac{3kc^{2}}{8\pi{G}\rho}}\cos{\alpha}$ then $$t+c=\Big(\frac{3kc^{2}}{8\pi{G}\rho}\Big)^{3/4}\frac{1}{\mp{i}\sqrt{kc^{2}}}\int^{\sqrt{8\pi{G}\rho/2kc^{2}}\arccos{a}}\sqrt{\cos{\alpha}}d\alpha=\Big(\frac{3kc^{2}}{8\pi{G}\rho}\Big)^{3/4}\frac{1}{\mp{i}\sqrt{kc^{2}}}E\Big(\frac{1}{2}\sqrt{8\pi{G}\rho/2kc^{2}}\arccos{a}\Big{|}2\Big)$$ Where $E(z, k^{2})$ is the incomplete elliptic integral of the second kind.