Borrar filtros
Borrar filtros

Problem solving a set of 5 non-linear ODE's

1 visualización (últimos 30 días)
tensorisation
tensorisation el 4 de En. de 2016
Comentada: tensorisation el 18 de En. de 2016
i have a problem when i try to solve a set of 5 odes. i tried using ode45 and the other ode's functions and it always gives me the following error:
Warning: Failure at t=8.822299e-22. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.009266e-36) at time t. > In ode45 (line 308)
here is my code:
c=2.988*10^10;
lambda=asin(0.9);
M_sun=1.988*10^33;
M=M_sun*10^9;
B=10^4;
G=6.67*10^-8;
r_H=G*M*(c^-2)*(1+cos(lambda));
Omega_H=(c^3/(2*G*M))*tan(lambda/2);
rho_rot=(Omega_H*B)/(4*pi*c);
rho_B=r_H;
Gamma_thr=10^6;
alpha=0.007297;
e=4.8032*10^-10;
C_1=(4/9)*alpha*Gamma_thr;
C_2=(e*Gamma_thr^4)/(4*pi*r_H^3*rho_rot);
m_e=9.10938*10^-28;
h_tilde_a=(Gamma_thr*m_e*c^2)/(4*pi*e*r_H^2*rho_rot);
function beta_a=beta_a(u_tilde_a)
beta_a=u_tilde_a.*(Gamma_thr^-2+u_tilde_a.^2).^(-1/2);
end
function Gamma_tilde_a=Gamma_tilde_a(u_tilde_a)
Gamma_tilde_a=(Gamma_thr^-2+u_tilde_a.^2).^(1/2);
end
function E_gamma_a=E_gamma_a(u_tilde_a)
E_gamma_a=(2/3)*(C_2/C_1)*Gamma_tilde_a(u_tilde_a).^3;
end
function alpha_tilde_a=alpha_tilde_a(u_tilde_a)
alpha_tilde_a=Gamma_tilde_a(u_tilde_a);
alpha_tilde_a(alpha_tilde_a<=1)=0;
alpha_tilde_a=alpha_tilde_a*C_1;
end
function q_tilde_a=q_tilde_a(n_tilde_lab_a,n_tilde_lab_plus,n_tilde_lab_minus,u_tilde_plus,u_tilde_minus)
q_tilde_a=(1./(2*n_tilde_lab_a)).*(n_tilde_lab_plus.*alpha_tilde_a(u_tilde_plus)+n_tilde_lab_minus.*alpha_tilde_a(u_tilde_minus));
end
function s_tilde_1_a=s_tilde_1_a(u_tilde_a)
s_tilde_1_a=C_2*(Gamma_tilde_a(u_tilde_a).^3).*u_tilde_a;
end
function q_tilde_1_a=q_tilde_1_a(n_tilde_lab_a,n_tilde_lab_plus,n_tilde_lab_minus,u_tilde_plus,u_tilde_minus)
q_tilde_1_a=(1./(2*n_tilde_lab_a)).*(n_tilde_lab_plus.*alpha_tilde_a(u_tilde_plus).*E_gamma_a(u_tilde_plus)+n_tilde_lab_minus.*alpha_tilde_a(u_tilde_minus).*E_gamma_a(u_tilde_minus));
end
function D_eqs=D_eqs(x,y)
D_eqs=zeros(5,1);
D_eqs(1)=y(2)-y(3)-1;
D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
end
s_tilde_initial=0;
s_tilde_final=100;
s_tilde_step=0.005;
beta_tilde_plus_0=-10^-3;
beta_tilde_minus_0=10^-3;
E_tilde_par_0=1;
n_tilde_lab_plus_0=abs(1/beta_tilde_plus_0);
n_tilde_lab_minus_0=abs(1/beta_tilde_minus_0);
u_tilde_plus_0=(1/Gamma_thr)*beta_tilde_plus_0*(1-beta_tilde_plus_0^2)^(-1/2);
u_tilde_minus_0=(1/Gamma_thr)*beta_tilde_minus_0*(1-beta_tilde_minus_0^2)^(-1/2);
[s_tilde,y_solved]=ode45(@D_eqs,[s_tilde_initial:s_tilde_step:s_tilde_final],[E_tilde_par_0 , n_tilde_lab_plus_0 , n_tilde_lab_minus_0 , u_tilde_plus_0 , u_tilde_minus_0]);
im at a loss here and clueless as to what is wrong here. i would very appreciate your help in the matter.
  22 comentarios
Mostrar 20 comentarios más antiguos
tensorisation
tensorisation el 11 de En. de 2016
Editada: tensorisation el 11 de En. de 2016
alpha is actually non-zero after the event, and zero before the event. but i also set a condition to check for the event inside D_eqs cuz untill the event has happened in each zone x>0 and x<0 2 of the 5 y functions are meaningless and i set them to zero there:
function D_eqs=D_eqs(x,y)
D_eqs=zeros(5,1);
D_eqs(1)=y(2)-y(3)-1;
if x>0
if Gamma_tilde_a(y(5))<=1
y(2)=0;
D_eqs(2)=0;
D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
y(4)=0;
D_eqs(4)=0;
D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
end
if Gamma_tilde_a(y(5))>1
D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
end
end
if x<0
if Gamma_tilde_a(y(4))<=1
D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
y(3)=0;
D_eqs(3)=0;
D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
y(5)=0;
D_eqs(5)=0;
end
if Gamma_tilde_a(y(4))>1
D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
end
end
end
i finally got it working but im unsure yet if the results make sense and there arent more problems to settle. i also wonder what if i have a case where the event is repeated multiple times: i start with f(y)<1, then the event first happen when f(y)=1, then f(y) goes up and eventually down and once again we have f(y)=1 and so on for a number of times? do i need to explicitly handle each event point, or is there an efficient way for dealing with this?
tensorisation
tensorisation el 18 de En. de 2016
now that i got it working, i ran into a strange error. i tried changing one of my parameters that appear in the equations from 1 to some higher value (say 2 for example) and i get the following error:
Error using horzcat Dimensions of matrices being concatenated are not consistent.
Error in ode15s (line 821) ieout = [ieout, ie];
Error in steady_model_cr (line 149)
im not sure what this error means, and what is the cause of it, since for exactly the same code i get no such error when one of my parameters is slightly lower.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (1)

John D'Errico
John D'Errico el 6 de En. de 2016
I won't try to run your code, nor will I try to figure out if there is a problem in the code, since I can't know if your code accurately reflects a realistic model. That is your job.
However, the error message indicates a classic problem - that the system is a stiff one. As importantly, the wide range of coefficients suggests this is HIGHLY probable to be a stiff system. Even the little I chose to read of your code suggests that stiff system is the VERY first thing I would consider.
Finally, I have no idea which tools you tried. Did you actually use one of the stiff solvers? If not, why not? Stiffness is the very first thing I would consider here, and ODE45 probably the last tool I'd throw at it, not the first. Note that the names of the stiff solvers in MATLAB all end in s.
  1 comentario
tensorisation
tensorisation el 6 de En. de 2016
Editada: tensorisation el 6 de En. de 2016
like i said in one of my comments, i tried using the other ode functions like ode15s and ode23s and i still get the same error.
as i've said in another comment, i tried running the code for a much more simple version of the equations, where the equations are much simpler and i also have only 1 constant (theres 2 constants there, that are exactly equal to each other):
h_tilde_a=Gamma_thr;
function Gamma_tilde_a=Gamma_tilde_a(u_tilde_a)
Gamma_tilde_a=(Gamma_thr^-2+u_tilde_a.^2).^(1/2);
end
function D_eqs=D_eqs(x,y)
D_eqs=zeros(5,1);
D_eqs(1)=y(2)-y(3)-1;
D_eqs(2)=( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1) );
D_eqs(3)=( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1) );
D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1) );
D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1) );
end
and i still get the same problem. so what can i possibly do to find out how to fix this problem?

Iniciar sesión para comentar.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by