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

조회 수: 7 (최근 30일)
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
Sam Chak
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;

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

채택된 답변

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

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by