ODE with a cubic root, ODE45 can't solve

Im trying to solve the following ODE, with has a peculair cubic-root as a nonlinearity (p). If this nonlinearity is uneven natural number, this is solved within a second.

But for the cubic root case, the solver doesn't seem to able to solve this. I've added a progress bar ( text-progress-bar-for-odes ) but it seems to be stuck at 0 %. Maybe the infinite derivative of x^1/3 at the origin is the problem.

Does anybody know how to get this simulation working? Thanks.

m = 1;
k = 1;
c = 0.000;
omega0 = sqrt(k/m);
x0_dot = 0.1;

Tl = 20;
%%NES

mna =0.02;
cna = 0.002;
kna = 0.4444;
p = 1/3;
 f1 = @(t,y)[y(3);y(4);...
    -(k*y(1) + c*y(3)+ kna*(y(1)-y(2))^p+cna*(y(3)-y(4)))/m;...
    -(kna*(y(2)-y(1))^p + cna*(y(4)-y(3)))/mna;...
    ];

Prec = 1e-14;

options = odeset('RelTol',Prec,'AbsTol',[Prec Prec Prec Prec],'OutputFcn',@odetpbar);

[T1,Y1] = ode45(f1,[0 Tl],[0 0 x0_dot 0],options);

Answer 1

[Torsten]

Instead of x^p, write sign(x)*(abs(x))^p.

Best wishes

Torsten.

[Steven Lord]

Or use nthroot. As an example, (-8)^(1/3) returns a complex value:

>> (-8)^(1/3)
ans =
   1.0000 + 1.7321i

There are three cube roots of -8, and one of them is real.

>> roots([1 0 0 8])
ans =
  -2.0000 + 0.0000i
   1.0000 + 1.7321i
   1.0000 - 1.7321i

The nthroot function returns the real cube root.

>> nthroot(-8, 3)
ans =
    -2

[dikkemulle]

Thanks! I should always check what something returns in Matlab before putting this in the ODE function!

Thanks,

Kevin