how to add a perturbation while solving ode using ode45? can i use the analytical approximation of Heaviside function for it if yes how?

조회 수: 2(최근 30일)
SOUVIK DARIPA 2023년 3월 1일
답변: Sam Chak 2023년 3월 3일
function dydt = rsf_ode4eq(t,y,d,par,v_lp,t_step,law)
dydt = zeros(4,1);
aa = par(1); bb = par(2); dc = par(3); sigma = par(4); stiff = par(5);
omega = y(1)*y(2)/dc;
dydt(1) = 1.d0 - omega;
%tau=1./(1+exp(-2*m*(t-d)));
% taudot =0.01*B_s*(2*m*exp(-2*m*(t-d))./((1+exp(-2*m*(t-d))).^2));
% if t <= t_step
% v_loc = v_lp;
% elseif t<=t_step+d
% v_loc=1.00001*v_lp;
% else
% v_loc = 1.00001*v_lp;
% end
if t <= t_step
v_loc = v_lp;
elseif t <= t_step+d
v_loc=1.000001*v_lp;
else
v_loc = 1.00000*v_lp;
end
% taudot_step = sirac(t,d,m,B_s);
% taudot_step=10e-8*B_s*(r-exp(-(t-d)/p));
% taudot_step=10-8*B_s*dirac(t-d);
% x=dirac(t-d);
tau_dot = stiff*(v_loc- y(2))+B_s*10e-6*x;
if strcmp(law,'A')
dydt(1) = 1.d0 - omega;
elseif strcmp(law,'S')
dydt(1) = -omega*log(omega);
end
dydt(3) = tau_dot;
end
solve:
clear
clc
close all
a = 0.021;
b = 0.025;
Dc = 100; % in microns
sigma = 1e1; % in MPa
mu_0 = 0.6;
m=10-9;
r=1;
p=1e-3;
Gmod = 3e4; % in MPa
b_s=mu_0*sigma;
v_shear = 3e9; % in microns/sec
Kc = sigma*(b-a)/Dc; % in MPa/microns
v_init = 1e-2; % in microns/s
theta_init = Dc/v_init; % in seconds
tspan = [1e-5 3.091924568069068e+07]; % in seconds
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.9*Im_stiff_osc; % Stiffness
tau_init=stiff*v_init;
B_s=mu_0*sigma;
law = 'A';
d=[2.5e7, 2.55e7, 2.6e7, 2.8e7 ,2.9e7 ,3e7];
for i = 1:2%length(d)
dd=d(i);
[t,y] = ode45(@(t,y)rsf_ode4eq(t,y,dd,par,v_init,t_step,law),...
tspan,[theta_init;v_init;tau_init;0],solopt);
theta=y(:,1);
v=y(:,2);
semilogy(t,v,'b','linewidth',0.75);
xlabel('Time[s]')
ylabel('Velocity[m/s]')
hold on
grid on
title('velocity vs time for perturbed system')
tau= sigma*(mu_0 + a*log(y(:,2)/v_init) + b*log(y(:,1)/theta_init));
hold on
% v_new(:,i)=v;
% t_new(:,i)=t;
% [pks,locs] = findpeaks(y(:,2),t);
% loc_new(:,i)=locs;
yyaxis right
semilogy(t,tau,'g','linewidth',0.75)
hold on
ylabel('\tau')
xlabel('time[s]')
title('\tau vs time[s] for perturbed system')
%
end
댓글 수: 3표시 이전 댓글 수: 2숨기기 이전 댓글 수: 2
SOUVIK DARIPA 2023년 3월 1일
to taudot @Sam Chak

댓글을 달려면 로그인하십시오.

채택된 답변

Sam Chak 2023년 3월 3일
If the perturbation looks like a step disturbance, then you can try using the Logictic function instead, which is a continuously differentiable. In physical systems, there will be a slight delay between the swiching states. Thus, I think it is reasonable to use that. Three parameters describe how the graph looks like:
x = linspace(-1, 2, 30001);
a = 3; % the supremum of the function
slope = 1e3; % steepness of the curve (adjust as you wish)
xsp = 0.5; % the point at which the switching occurs
approx_step = a./(1 + exp(- slope*(x - xsp))); % a logistic function
plot(x, approx_step, 'linewidth', 1.5), grid on
xlabel('x'), ylim([-1 4])
Then, try adding perturbation like this:
dydt(1) = 1.d0 - omega;
dydt(3) = tau_dot + approx_step;

댓글을 달려면 로그인하십시오.

범주

Find more on Mathematics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by