ODE15s: segmentation violation and imaginary parts

Hi. I am solving

If I set P0=2100 I get a result, while if I set P0=2000: the first time I run the code, MATLAB returns "Segmentation violation"; if I try a second time running the code, MATLAB returns "Warning: Imaginary parts of complex X and/or Y arguments ignored." Can anyone help me out? Thank you.

function out = convezione_20(t,y,R0)
S = 0.073;
rho = 998;    
mi = 1.005e-3;
P0 = 2000; %<-------------------------------------------
L = 2454000;
Tinf = 293;   
rho_g = 1/57.79;
lambda = 0.5982;
cp = 4184;    
Pv_Tinf = 2329.6;
%
out=[y(2)
    -1.5 * y(2)^2 / y(1) + 1 / rho / y(1) * (Pv_Tinf - P0 +...
    (L * rho_g)^2 / Tinf * y(2) / y(3) -...
    4 * mi * y(2) / y(1) - 2 * S / y(1))
    -y(3) + lambda * 0.5 / y(1) * (2 + (-6 * rho * cp * y(2) / pi / y(3))^(1/3)...
    - 12 * rho * cp * y(2) / pi / y(3)) ];

run file

R0 = 510e-6;
y0 = [R0; 0; 1172.94];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M);%,'RelTol',1e-10,'AbsTol',[1e-10 1e-10 1e-10]);
[t,y] = ode15s(@(t,y) convezione_20(t,y,R0),[0 40e-6],y0,options);
t = t * 1000000;
y = y * 1000000;
[t,y(:,1)];
plot(t,y(:,1))

Edit: this is what MATLAB says:

(...)

Register State:
  EAX = 00000008  EBX = 00000000
  ECX = 02723d90  EDX = 00000000
  ESI = 00cdb99c  EDI = 00cdb988
  EBP = 00000000  ESP = 00cdb750
  EIP = 00000000  FLG = 00210246

(...)

Error in ==> convezione_20 at 13
out=[y(2)

Error in ==> run_convezione_20>@(t,y) convezione_20(t,y,R0) at 5
[t,y] = ode15s(@(t,y) convezione_20(t,y,R0),[0 40e-6],y0,options);

Error in ==> funfun\private\odenumjac at 131
      Fdel(:,j) = feval(F,Fargs{1:diffvar-1},ydel(:,j),Fargs{diffvar+1:end});

Error in ==> funfun\private\daeic12 at 38
  [DfDy,Joptions.fac,nF] = odenumjac(fun, {t0,y,args{:}}, f, Joptions);  

Error in ==> ode15s at 374
    [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in ==> run_convezione_20 at 5
[t,y] = ode15s(@(t,y) convezione_20(t,y,R0),[0 40e-6],y0,options);

[Torsten]

Check whether -6*rho*cp*Rdot/(pi*h) becomes negative.

As a consequence, (...)^1/3 gets complex.

Best wishes

Torsten.

[Marco Sammito]

Why? If n<0, n^(1/3) is not complex.

[Torsten]

Try it.

http://www.wolframalpha.com/input/?i=(-1)%5E(1%2F3)

Best wishes

Torsten.

[Marco Sammito]

I tried using as exponent 1 insted of 1/3, but now the warning becomes "Warning: Failure at t=2.874099e-007. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.021085e-021) at time t."

[Torsten]

Does your initial value for h (h0) together with the initial values for R and Rdot (R0 and Rdot0) satisfy the algebraic equation, i.e.

-h0 + lambda * 0.5 / R0 * (2 + (-6 * rho * cp * Rdot0 / pi / h0)^(1/3) - 12 * rho * cp * Rdot0 / pi / h0) = 0

?

If not, first try to adjust h0 such that this condition holds.

Best wishes

Torsten.

[Jan]

@Marco: Here are 2 problems:

1. The segementation violation
2. The mathematical problem

It is not worth to examine the 2nd one before the 1st is solved, because a segmentation violation means, that either the RAM, the processor, Matlab itself or a MEX function works incorrectly. You cannot solve the mathematical problem, when the computing system is not working reliably.

See [EDITED in my answer]

Answer 1

[Jan]

A segmentation violation is severe. You cannot trust any computations afterwards. I cannot imagine that the posted code produces such an error. Do you use any user-defined Mex-functions?

[EDITED] The segmentation violation must be solved at first:

• Check again in any user defined C-Mex code is running:

[mFun, mexFun] = inmem

after your code ran.

• Run exhaustive unit-tests, if you have some.
• Run examples from the documentation with ode15s.
• Check the hardware: CPU, memory, ...
• Re-install Matlab

When the system is running reliably, caring about the mathematical problem is the next step. The input parameters of the mathematical problem might trigger the segmentation fault, but they are not the reason.

[Marco Sammito]

I don't have any Mex-functions. What do you mean with "You cannot trust any computations afterwards"? Even if, say, I turn the computer off and fire it up again?

[Jan]

A segmentation vialoation means, that some code has written to a wrong memory address. This could cause any arbitrary behavior. Restarting Matlab is obligatory (and as far as I know this happens automatically). But if the code runs the next time, this might be luck only: Code, which creates memory violations is not trustworthy.

Do you get any further information? "MATLAB returns 'Segmentation violation'" is a lean description of a serios error.

[Marco Sammito]

I edited the original post.