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 (최근 30일)
이전 댓글 표시
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
Torsten
2023년 4월 14일
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.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!