If there any possible solution (analytical/numerical) of a this nonlinear third order differential equation?

Question

The equation I'm trying to solve is the following $$m_3\ddot{y}^3 + m_1\ddot{y} + k_1 y + k_3y^3 = f_1(t) $$

All the parameters are real $m_1,m_3,k_1,k_3\in \mathbb{R}$ and constant. $f_1(t)\in \mathbb{R}$ is a function of $t$ as is $y(t)$. Obviously if $m_3=0$ then the solution can be found because the equation takes the form

$$\ddot{y} = a y + by^3 + c,$$

where $a=-\dfrac{k_1}{m_1}$, $b = -\dfrac{k_3}{m_1}$, $c = \dfrac{f_1}{m_1}$, which is the classical form of a nonlinear ODE as $$ \ddot{y} = F(y,t).$$

EDIT:

Please find below my implementation in MATLAB of DAE with NO success so far.

First I define the DAE function as

function out = funTest(t,y,yp)

% coefficients
m1 = 9.8845614144e-10;
m3 = 2.00030511601069e-12;

f1 = 1.66666666666667e-05;

k1 = 0.214188819381172;
k3 = 8.6610733212416;

omega = 2*pi*(-82.454*t + 10030);

out = [yp(1) - y(2)
    yp(2) - y(3)
    m3*y(3)^3 + m1*y(3) + k1*y(1) + k3*y(1)^3 - f1*sin(omega*t)];

Below the main file

tspan = [60 120];
M = [1 0 0; 0 1 0; 0 0 0]; % constant mass matrix
options = odeset('RelTol',1e-4,'Jacobian',{[],M});

y0 = [0; 0; 0];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@funTest,0,y0,[1 1 0],yp0,[],options);

[t,y] = ode15i(@funTest,tspan,y0,yp0,options);


plot(t,y(:,1));
ylabel('y');

I tried changing both options as

options = odeset('RelTol',1e-4,'AbsTol',[1e-10 1e-6 1e-6], ...
   'Jacobian',{[],M});

but nothing happens

Answer 1

Let's discuss the scales of the example, where $m_1,m_3,k_1,k_3,f_1$ have scales $10^{-9}, 10^{-12}, 1, 1, 10^{-5}$ respectively. The "easy" adjustments are changing the time scale and the scale of $y$. Set $y(t)=10^{-m}u(10^{n}t)$, then $\ddot y=10^{2n-m}\ddot u(10^{n}t)$ and the transformed equation is $$ 10^{6n-3m}m_3\ddot u(s)^3+10^{2n-m}m_1\ddot u(s)+10^{-3m}k_3u(s)^3+10^{-m}k_1u(s)=f(t) $$ Now to balance the size of the coefficients, multiply with $10^{m}$ to get $$ 10^{6n-2m}m_3\ddot u(s)^3+10^{2n}m_1\ddot u(s)+10^{-2m}k_3u(s)^3+k_1u(s)=10^{m}f_1\sin(\theta(t)), $$ and then identify as the active constraints $m\sim 5$ for the right side and $6n\lesssim 12+2m$ and $2n\lesssim 9$ for the second derivatives. One variant is $m=6$, $n=4$. Then the transformed constants $m_1,m_3,k_1,k_3,f_1$ have the sizes $10^{-1}, 1, 1, 10^{-12}, 10$. The dominant terms now have sizes that are compatible with the internal heuristics and default tolerances of the solvers. One can now translate this rescaling back into sensible absolute tolerances for the problem in the original coefficients. For a goal of 4 valid digits, $y$ would have AbsTol $10^{-m-4}=10^{-10}$, $\dot y$ of $10^{n-m-4}=10^{-6}$ and the equation of $10^{-5-4}=10^{-9}$.

Numerical method as ODE solving the cubic in each step

Use the balanced scaling computed above. To solve the cubic $m_3a^3+m_1a=c$ a first approximation is in general $a=(c/m_3)^{1/3}$. For very small $c$ it is $a=c/m_1$. Both behaviors are combined in the formula $$ a = \frac{c}{(m_1^3+m_3c^2)^{1/3}},~~\text{ or }~~a = \frac{c}{\max(m_1,(m_3c^2)^{1/3})}, $$
which additionally is unproblematic for negative $c$. This gives a very good initial value for the numerical solver fsolve.

from numpy import sin,pi
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

# coefficients
m1 = 9.8845614144e-10*1e8;       # *10^(2n), m=6, n=4
m3 = 2.00030511601069e-12*1e12;  # *10^(6n-2m)

f1 = 1.66666666666667e-05*1e6;   # *10^m

k1 = 0.214188819381172;          # unchanged
k3 = 8.6610733212416*1e-12;      # *10^(-2m)

s0 = 60*1e4; sf = 120*1e4;

# angle function of the forcing
def theta(s): 
    t = 1e-4*s; # s = 10^n*t
    omega = 2*pi*(-82.454*t + 10030); 
    return omega*t;

# system function
def ODEsys(s,u):
    c = f1*sin(theta(s)) - k1*u[0]-k3*u[0]**3;
    du = u[1];
    ddu = fsolve(lambda z: m1*z+m3*z**3-c, c/(m1**3+m3*c**2)**(1./3))
    return [du, ddu]

# compute a small initial segment of the solution
sol = solve_ivp(ODEsys, [s0,s0+1000], [0,0], method='Radau', max_step=1)
# plot in the original scale, t = 10^(-n)*s, y(t)=10^(-m)*u(s)
plt.plot(sol.t*1e-4, sol.y[0]*1e-6);
plt.grid(); plt.show()

which gives the plot

For the slightly larger segment to $t=61$ this method gives the plot

plots for larger time scales like to $t=120$ do not appear sensible.

Answer 2

The problem is to avoid the potential complex solutions for the cubic. Using a symbolic processor we can follow with the cubic real root (always exists). Using MATHEMATICA as the symbolic solver we submit

sols = Solve[m3 y''[t]^3 + m1 y''[t] + k1 y[t] + k3 y[t]^3 == f[t],y''[t]]

obtaining as real root

$$ y''(t) \to \frac{\sqrt[3]{2} \left(9 m_3^2 \left(f(t)-y(t) \left(k_3 y(t)^2+k_1\right)\right)+\sqrt{3} \sqrt{m_3^3 \left(27 m_3 \left(-f(t)+k_3 y(t)^3+k_1 y(t)\right){}^2+4 m_1^3\right)}\right){}^{2/3}-2 \sqrt[3]{3} m_1 m_3}{6^{2/3} m_3 \sqrt[3]{9 m_3^2 \left(f(t)-y(t) \left(k_3 y(t)^2+k_1\right)\right)+\sqrt{3} \sqrt{m_3^3 \left(27 m_3 \left(-f(t)+k_3 y(t)^3+k_1 y(t)\right){}^2+4 m_1^3\right)}}} $$

and after that we solve numerically the DE

m1 = 9.8845614144 10^-10;
m3 = 2.00030511601069 10^-12;
f1 = 1.66666666666667 10^-05;
k1 = 0.214188819381172;
k3 = 8.6610733212416;
omega = 2*Pi*(-82.454*t + 10030);
f[t_]:= f1*Sin[omega*t]
tmax = 1;
equ = y''[t] == (y''[t]/.sols[[1]]);
soly = NDSolve[{equ, y[0] == 0.1, y'[0] == 0.1},y,{t,0,tmax}][[1]];
Plot[Evaluate[y[t] /. soly], {t, 0, tmax}]