Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab ODE15s error message

  1. Apr 21, 2005 #1

    Clausius2

    User Avatar
    Science Advisor
    Gold Member

    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 MATLAB:eek:de15s: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
     
  2. jcsd
  3. Apr 21, 2005 #2
    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 occured 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.
     
  4. Apr 21, 2005 #3

    Clausius2

    User Avatar
    Science Advisor
    Gold Member

    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 [tex]F(\eta)[/tex] and [tex]Y(\eta)[/tex] which represents the self similar stream function and mass fraction respectively. The system obtained is:

    [tex] (F'Y)''+\frac{F}{2}(F'Y)'=0 [/tex] (1)

    [tex] \Big(\frac{Y'}{Y}\Big)'+\frac{F}{2}Y'=0 [/tex] (2)


    Boundary conditions are:

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

    2) [tex]\eta\rightarrow-\infty[/tex] then [tex] (F'Y)=Y=1[/tex] and [tex]F=\eta[/tex]

    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.
     
    Last edited: Apr 21, 2005
  5. Apr 21, 2005 #4
    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.
     
  6. Apr 21, 2005 #5
    What is f defined as in the failure line? My understanding is that f is a seperate 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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Matlab ODE15s error message
Loading...