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

3 views (last 30 days)
tensorisation
tensorisation on 4 Jan 2016
Commented: tensorisation on 18 Jan 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 Comments
Show 20 older comments
tensorisation
tensorisation on 11 Jan 2016
Edited: tensorisation on 11 Jan 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 on 18 Jan 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.

Sign in to comment.

Sign in to answer this question.

Answers (1)

John D'Errico
John D'Errico on 6 Jan 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 Comment
tensorisation
tensorisation on 6 Jan 2016
Edited: tensorisation on 6 Jan 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?

Sign in to comment.

Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!