# Matlab: ode45 loop

I have to solve this equation

$$\ddot{x}x = -\frac{3}{2}\dot{x}^2 + \frac{\dot{x}}{h} + \sin(t)$$

where $h$ is defined by

$$h = \left(\frac{\dot{x}}{h}\right)^{1/3} + \frac{\dot{x}}{h}$$

My idea is to initialize $h=1$, solve the first equation, calculate a new value for h:

$$h = \frac{\dot{x}}{x\ddot{x} + (3/2)\dot{x}^2 - \sin(t)}$$

and do this many times because it guarantees that h has converged. In the following code I do not know how to express $\ddot{x}$, in fact MATLAB says that xp ($\dot{x}$) is undefined.

Plesset.m file

function xp = Plesset(t, x)
xp = zeros(2, 1);
xp(1) = x(2);
h = 1;
h = (x(2)/h)^(1/3) + x(2)/h;
xp(2) = 1/x(1)*(-1.5*(x(2))^2 + x(2)/h+sin(t));


run.m file

for i = 1:100
[t, x] = ode45('Plesset', [0,5], [0,0]);
h = x(2)/(x(1)*xp(2) + 3/2*(x(2))^2 - sin(t));
end
[t, x(:,1)]
plot(t, x(:,1))


Let's first take your equation for $h$, which is defined as:

$$h = \left(\frac{\dot{x}}{h}\right)^{1/3} + \frac{\dot{x}}{h}$$

If we define $u=\frac{\dot{x}}{h}$, we can transform the equation to:

$$0 = u^{2} + u^{4/3} - \dot{x}$$

As you can notice, this equation can be seen as close to a parabola with respect to $u$, where the difference is the $u^{4/3}$ grows a little faster than linearly. To find the approximate parabolic form of this equation, we can use the Taylor Series of $u^{4/3}$ about $u=1$, which gives:

$$u^{4/3} \approx 1 + \frac{4}{3}(u-1) + \frac{2}{9}(u-1)^2$$

Substitute this quantity into the transformed equation and simplify and you'll get:

$$11u^2 + 8u - (9\dot{x}+1) \approx 0$$

Solving this equation for $u$, we get:

$$u \approx \frac{1}{11} \left(-4 \pm \sqrt{16 + 11(9\dot{x}+1)}\right)$$

This approximate solution for $u$ should be fairly close to the actual solution, though you could refine it using one or two Newton iterations. The recursive equation for a single Newton Iteration is:

$$u = u - \frac{u^2 + u^{4/3} - \dot{x}}{2u + \frac{4}{3}u^{1/3}}$$

Note that I use the most positive root the approximation produces. In your dynamics, you could then actually use $u$ directly using:

$$\ddot{x} = \frac{1}{x}\left(\sin(t) + u - \frac{3}{2}\dot{x}^2\right)$$

Using the approximation for $u$, without using any Newton iterations, and the initial conditions $[x,\dot{x}]^{T} = [1,0]^{T}$ produces the following results using ode45(...):

• Why ode45? And the analytic manipulations are to get a better guess? – Pu Zhang Nov 9 '16 at 3:03
• @PuZhang I used ode45 because it will work fine and the OP used it. The analytic manipulations lead to a very good guess that closely approximates the correct solution. So with only a couple Newton iterations you can achieve great precision. – spektr Nov 9 '16 at 3:08
• @choward With the initial conditions [x, \dot{x}]=[0,1] Matlab returns NaN. Instead, with the initial conditions [x, \dot{x}]=[1,0.001] I got the result. And why did you use the Taylor series about u=1? – Marco Nov 9 '16 at 8:49
• @Marco You get NaN because of the $\frac{1}{x}$ term in the equation for $\ddot{x}$. You can't start at $x = 0$ based on the differential equation. And I did the Taylor Series about one so the first and second derivative terms would evaluate to 1, so no special reasoning. I knew the series would be close either way since that term is superlinear. – spektr Nov 9 '16 at 13:23

You should adopt a different approach here:

1. Express your second order equation as a system of two first order equations in

$u(t)=\left[ \begin{array}{c} x(t) \\ \dot{x}(t) \end{array} \right]$

1. Your system of equatrions will be "implicit." That is, you can write it as

$f(\dot{u}(t), u(t), t)=0$

MATLAB has the ability to work directly with implicit systems of differential equations in this form. See the ode15i solver.

• But I have already expressed it as a system of two first order equations: xp(1)=x(2); and xp(2)=1/x(1)*(-1.5*(x(2))^2+x(2)/h+sin(t)); – Marco Nov 8 '16 at 17:35