How do I solve this Differential Algebraic Equation?

Question

I am trying to model a physical system which is goverened by the equation:

$$ (r(t))^2 \cdot \frac{d^2 r(t)}{dt^2} = c q(t)$$

I am a junior electrical engineering student, so I have had courses in linear algebra, ordinary differential equations, and LTI system theory + Fourier/Laplace/DFT/Z transforms.

Based on what I can tell about this equation, it is a Differential Algebraic Equation, which I am unfamiliar with how to solve. I would like to solve for $r(t)$ given an arbitrary function of $q(t)$ - say a sinusoidal, exponential, or polynomial function.

From what I've read online, an analytic solution may not be possible, so numerical solutions are fine too.

Also, all I can find online is how to solve a SYSTEM of DAEs, not simply a singular DAE equation.

What I would like to do is to be able to predict the output $r(t)$ given some input signal $q(t)$

Answer 1

As Yves mentioned your problem is an ordinary differential equation, however it is implicit (the highest derivative is not isolated). If you are expecting the solution cross the zero line the reformulation into an explicit system by diving by $r(t)^2$ is inappropiate. I will focus my answer on the implicit formulation.

Please note that $\dfrac{q(t)}{r^2(t)}$ is not Lipschitz-continous (has unbounded derivative) meaning it has likely no solution for all points in time/value space. Solvers might fail to solve it over the complete time interval you want it to be solved. You can solve implicit ODEs numerically using an DAE solver. The idea is that DAE solver tries to find an du such that the residual function returns 0 and then step forward using that du. Here is how you would do it in a Julia command line:

]add DifferentialEquations #installs a kitchen sink of tools
]add Plots # if you need plots

function residual(res,du,u,p,t)
   res[1] = du[1] - u[2]
   res[2] = du[2]*u[1]*u[1] - c*q(t)
end

c = 3 # when changed just rerun solve
q(t) = sin(t) # when changed just rerun solve

tspan = (0.0,1.0) #solve for virtual 1.0 time units
u0 = [1.,2.] # initial conditions for r(t) and s(t)

prob = DAEProblem(DAEFunction(residual),u0,u0,tspan) # If you change u0, or tspan rerun this too

sol = solve(prob, DABDF2())

sol(0.8) #interpolated solution at time 0.8
Array(sol) #solution without interpolation

plot(sol) ## savefig(plot(sol), "name.png") to save in working directory, you can use a path too

Another reformulation would be using a mass matrix. The idea would be that you bring your equation into this form. $$M(\vec{x},t) \dfrac{d\vec{x}}{dt} = f(\vec{x},t)$$

However whether you have problems at zero crossings (where the $M^{-1}$ is illdefined) depends on the solver in question, i think. If you want to see code for that too, please add a comment saying so. I could also show how to solve the explicit formulation of Yves in Julia.