Problem solving a set of 5 non-linear ODE's
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
John D'Errico
le 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 commentaire
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!