how can I dissect a ode in 3 parts ! lets say I am starting from [0 t1] then stop the ode add a value to a parameter again start the ode[ t1 t2 ]then again change it back ?
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
sd()
%this is the code, lets say i dont add anything to the code then i should get back the original cycles
function sd
a = 0.021;
b = 0.025;
Dc = 100; % in microns
sigma = 1e1; % in MPa
mu_0 = 0.6;
Gmod = 3e4; % in MPa
v_shear = 3e9; % in microns/sec
rad = 0.5*Gmod/v_shear; % Radiation damping term
v_dyn = (a*sigma)/rad;
Kc = sigma*(b-a)/Dc; % in MPa/microns
v_init = 1e-2; % in microns/s
theta_init = Dc/v_init; % in seconds
tspan1=[0 2.5e7];
solopt = odeset('RelTol',1e-16);
t_step = 1e5; % in seconds
Im_stiff_osc = sigma*(sqrt(b)-sqrt(a))^2/Dc;
rat_oscillatory_sol = Im_stiff_osc/Kc;
stiff = 0.7*Im_stiff_osc; % Stiffness
B_s=mu_0*sigma;
law = 'A';
par = [a,b,Dc,sigma,stiff,rad];
taudot_step=[0 9*10e-9*B_s 2.6*10e-9*B_s];
[t1,y11] = ode45(@(t,y)ode5(t,y,par,v_init,t_step,law),...
tspan1,[theta_init;v_init],solopt);
theta_init=y11(end,1);
v_init=y11(end,2);
tspan2=tspan1+2.51e7;
[t2,y12]=ode45(@(t,y)ode5(t,y,par,v_init,t_step,law),...
tspan2,[theta_init;v_init],solopt);
figure()
semilogy(t1,y11(:,2),'r')
hold on
semilogy(t2,y12(:,2),'g')
ylabel ('$velocity$', 'Interpreter','latex')
xlabel ('time[s]')
hold off
end
function dydt = ode5(t,y,par,v_lp,t_step,law)
aa = par(1); bb = par(2); dc = par(3); sigma = par(4); stiff = par(5);
rad = par(6);
dydt = zeros(2,1);
omega = y(1)*y(2)/dc;
dydt(1) = 1.d0 - omega;
if t <= t_step
v_loc = v_lp;
elseif t > t_step
v_loc=1.000001*v_lp;
else
v_loc = 1.00000*v_lp;
end
tau_dot = stiff*(v_loc- y(2));
if strcmp(law,'A')
dydt(1) = 1.d0 - omega;
elseif strcmp(law,'S')
dydt(1) = -omega*log(omega);
end
dydt(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);
end
4 commentaires
Torsten
le 14 Avr 2023
Yes, of course you will get the same plot if you don't change anything in the equations. But there seem to be singularities in your solution function.
Réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!