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

Question

SOUVIK DARIPA 2023-3-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1921095-how-to-add-a-perturbation-while-solving-ode-using-ode45-can-i-use-the-analytical-approximation-of

回答： 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);
rad = par(6);B_s=par(7);d=par(8);m=par(9);r=par(10);p=par(11);
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(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);
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
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
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;
% xlsread('t_0.xlsx');
% r_numf=xlsread('t_0','C2:C52');
law = 'A';
d=[2.5e7, 2.55e7, 2.6e7, 2.8e7 ,2.9e7 ,3e7];
par = [a,b,Dc,sigma,stiff,rad,B_s,d,m,r,p];
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 个评论
显示 1更早的评论隐藏 1更早的评论

Sam Chak 2023-3-1

A total of 98 lines in the code without explanatory comments.

Where exactly do you want to add the perturbation?

dydt(1) = 1.d0 - omega;
dydt(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);
dydt(3) = tau_dot;

SOUVIK DARIPA 2023-3-1

to taudot @Sam Chak

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2023-3-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1921095-how-to-add-a-perturbation-while-solving-ode-using-ode45-can-i-use-the-analytical-approximation-of#answer_1184510

Hi @SOUVIK DARIPA

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(2) = ((tau_dot/sigma - bb*dydt(1)/y(1)))/(aa/y(2)+rad/sigma);

dydt(3) = tau_dot + approx_step;