Troubleshooting Matlab ODE15s Error Message: Understanding f and Options

[Clausius2]

I am employing the Ode15s subroutine of Matlab for solving some differential equation. But I have obtained the next error message:

"Warning: Failure at t=-6.718807e-001. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.387000e-015) at time t.
(Type "warning off MATLABde15s:IntegrationTolNotMet" to suppress this warning.)
> In C:\MATLAB6p5\toolbox\matlab\funfun\ode15s.m at line 743"

I don't know what does it mean, and how can it be solved. If I change the ode solver it remains happening the same.

Thanks

 

[Lonewolf]

What differential equation are you trying to solve? I think it's either going to be your function is not defined everywhere. For example, you've got a square root or a log and discretisation errors have occurred causing the input to become negative.
Your equation might also have a singularity somewhere, as in y' = y^2, y(0)=1, which is indeed singular at t=1.

 

[Clausius2]

What differential equation are you trying to solve? I think it's either going to be your function is not defined everywhere. For example, you've got a square root or a log and discretisation errors have occurred causing the input to become negative.
Your equation might also have a singularity somewhere, as in y' = y^2, y(0)=1, which is indeed singular at t=1.

Thanks for your reply. This is my problem:
In order to solve the near field description of a round jet, I want to work out the variables $$F(\eta)$$ and $$Y(\eta)$$ which represents the self similar stream function and mass fraction respectively. The system obtained is:

$$(F'Y)''+\frac{F}{2}(F'Y)'=0$$ (1)

$$\Big(\frac{Y'}{Y}\Big)'+\frac{F}{2}Y'=0$$ (2)


Boundary conditions are:

1) $$\eta\rightarrow\infty$$ then $$F=C_1+C_2e^{-0.5C_1\eta}(\eta+C_3+2/C_1)$$ and $$Y=\frac{1}{0.5C_1(\eta+C_3)}$$; ; ;

2) $$\eta\rightarrow-\infty$$ then $$(F'Y)=Y=1$$ and $$F=\eta$$

I have to work out the constants C1, C2 and C3 shooting from +infinite and yielding the boundary conditions at -infinite. My program is:
__________________________________________________________
function f=f(eta,y);

%y1=F; y2=F'/rho; y3=(F'/rho)'; y4=Y; y5=rho*Y'; rho=1/Y;

%Sistema equivalente:

%y1'=rho*y2;
%y2'=y3;
%y3'=-0.5*y1*y3;
%y4'=y5/rho;
%y5'=-0.5*Sc*y1*y5/rho;

f(1)=y(2)/y(4);
f(2)=y(3);
f(3)=-0.5*y(1)*y(3);
f(4)=y(5)*y(4);
f(5)=-0.5*y(1)*y(5)*y(4);
f=f';
__________________________________________________

clear; clc; clf; format long;

%INICIALIZACION DE CONSTANTES
%----------------------------------
C1=2.18; %F_inf
C2=1e-3; %A
C3=1e-3; %eta_inf

C=[C1 C2 C3];

%PARAMETROS DE INTEGRACION
%-------------------------------
maxerror=1e-4;
e=0; %relacion de pesos moleculares;
S=1; %#Schmidt
error=10*maxerror;
iter=0;
Z=0.01;

%PLANO COMPUTACIONAL
%-------------------------------
etamax=10;
etamin=-20;
etao=[etamax etamin];
h=0.1;
neta=(etamax-etamin)/h+1;
%etao=linspace(etamax,etamin,neta);

%INTEGRACION
%--------------------------

pause;
while abs(error)>maxerror
iter=iter+1;
Y0(1)=C(1)+C(2)*exp(-0.5*C(1)*etamax)*(etamax+C(3)+2/C(1)); %F
Y0(4)=(0.5*S*C(1)*(etamax+C(3)))^(-1); %Y
Y0(2)=C(2)*exp(-0.5*C(1)*etamax)*(etamax*(1-0.5*C(1))-0.5*C(1)*C(3)-1)*Y0(4); %(F'Y)
Y0(3)=-Y0(2)/(etamax+C(3))+C(2)*exp(-0.5*C(1)*etamax)*(1-0.5*C(1))*Y0(4)-0.5*C(1)*Y0(2); %(F'Y)'
Y0(5)=-1/(etamax+C(3)); %(Y'/Y)

%options = odeset('RelTol',1e-4,'AbsTol',1e-4);
options=[];

[eta,y]=ode15s('f',etao, Y0);<------- HERE'S THE FAILURE LINE
N=length(eta);
y=real(y);

R(1)=y(N,1)-etamin;
R(2)=y(N,2)-1;
R(3)=y(N,4)-1;
error=sum(abs(R));

pause;

for i=1:3
C1=C;
dC=Z*C(i);
C1(i)=C1(i)+dC;

Y01(1)=C1(1)+C1(2)*exp(-0.5*C1(1)*etamax)*(etamax+C1(3)+2/C1(1)); %F
Y01(4)=(0.5*S*C1(1)*(etamax+C1(3)))^(-1); %Y
Y01(2)=C1(2)*exp(-0.5*C1(1)*etamax)*(etamax*(1-0.5*C1(1))-0.5*C1(1)*C1(3)-1)*Y01(4); %(F'Y)
Y01(3)=-Y01(2)/(etamax+C1(3))+C1(2)*exp(-0.5*C1(1)*etamax)*(1-0.5*C1(1))*Y01(4)-0.5*C1(1)*Y01(2); %(F'Y)'
Y01(5)=-1/(etamax+C1(3)); %(Y'/Y)

options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[eta1,y1]=ode15s('f',etao, Y01);
N1=length(eta1);
y1=real(y1);

R1(1)=y1(N1,1)-etamin;
R1(2)=y1(N1,2)-1;
R1(3)=y1(N1,4)-1;

A(1,i)=(R1(1)-R(1))/dC;
A(2,i)=(R1(2)-R(2))/dC;
A(3,i)=(R1(3)-R(3))/dC;

end

delta=inv(A)*R';
lambda=-1; %damping factor
delta=(lambda*delta);
C=C+delta';
disp([C iter error]);
figure(1); subplot(2,2,1); plot(y(:,1),eta); xlabel('F');
subplot(2,2,2); plot(y(:,2),eta); xlabel('F´/\rho');
subplot(2,2,3); plot(y(:,4),eta); xlabel('Y');
pause;
end

___________________________________________________-

Don't worry about spanish words. Look at the failure line I have bolded. The error points that line. The code is unable to integrate the problem at the first iteration. I don't think the variables are discontinous. I don't know the problem.

 

[Lonewolf]

I can't see any obvious errors. I'd try altering f and see if you get the same errors. Got me stumped to be honest.

 

[physicsCU]

What is f defined as in the failure line? My understanding is that f is a separate function file, or needs to be.

Also, the line above it with options could be a trouble spot.

Those are the two spots I would look at.

Oh, try putting an @ in front of f and removing the ' marks. Just a syntax thing, but it will help.

Those are the three things I would look at first. If none of those work, I am stumped too.