Lonewolf said:
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.
de15s:IntegrationTolNotMet" to suppress this warning.)