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));
[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(...):

enter image description here

| cite | improve this answer | |
  • $\begingroup$ Why ode45? And the analytic manipulations are to get a better guess? $\endgroup$ – Pu Zhang Nov 9 '16 at 3:03
  • $\begingroup$ @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. $\endgroup$ – spektr Nov 9 '16 at 3:08
  • $\begingroup$ @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? $\endgroup$ – Marco Nov 9 '16 at 8:49
  • $\begingroup$ @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. $\endgroup$ – 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.

| cite | improve this answer | |
  • $\begingroup$ 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)); $\endgroup$ – Marco Nov 8 '16 at 17:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.