dirac() and ode45 issue

10 views (last 30 days)
Ismael Ochoa
Ismael Ochoa on 22 Oct 2021
Answered: Star Strider on 22 Oct 2021
I am trying to solve an ODE that includes the direct function, but I don't see the efect of it whe solving the ODE. I have reviewed some previous anwswers to how to use direc() function and ODE, I have followed a recomended steps: run the ODE till the point when the dirac() is triggered, get the final values of the variables and run again the ode45 using those values as initial condition ( see below my code.). However the first ODE output is "inf" and I don't see any output from the second ODE.
clear;
close;
tspan1 = 0:0.00000001:0.01;
x0=[0 0];
%first ODE within 0 to 0.01, 0.01 DIRAC() is triggered.
[t1,x1] = ode45(@cirt,tspan1,x0);
xdirac = x1(length(t1),1); %initial condition for the next ODE
tdirac = t1(length(t1),1); %initial condition for the next ODE
% next ODE
tspan = 0.01:0.000001:0.02;
[ttt,x] = ode45(@cirt,tspan,[tdirac xdirac]);
plot(ttt,x(:,1))
grid on;
function dxdt=cirt(t,x)
r1 = 0.1;
r2 = 0.1;
R=300;
L=0.047;
C=0.33E-6;
Vs = heaviside(t-0.01);
b = 1/(R*C);
a2 = L*(1+r2/R);
a1 = L/(R*C) + r1*(1+r2/R)+r2;
a0 = r1/(R*C) + 1/C;
dx1dt = x(2);
dx2dt = (-a1 / a2 )* x(2) + (-a0/a2)*x(1) + b*Vs/a2 - (1+r2/R) * dirac(t-0.01);
dxdt = [dx1dt ;dx2dt];
end

Answers (1)

Star Strider
Star Strider on 22 Oct 2021
I doubt that the dirac call is actually necessary.
The code works correctly without it —
tspan1 = 0:0.00000001:0.01;
x0=[0 0];
%first ODE within 0 to 0.01, 0.01 DIRAC() is triggered.
[t1,x1] = ode45(@cirt,tspan1,x0);
xdirac = x1(length(t1),1); %initial condition for the next ODE
tdirac = t1(length(t1),1); %initial condition for the next ODE
% next ODE
tspan = 0.01:0.000001:0.02;
[ttt,x] = ode45(@cirt,tspan,[tdirac xdirac]);
plot(ttt,x(:,1))
grid on;
function dxdt=cirt(t,x)
r1 = 0.1;
r2 = 0.1;
R=300;
L=0.047;
C=0.33E-6;
Vs = heaviside(t-0.01);
b = 1/(R*C);
a2 = L*(1+r2/R);
a1 = L/(R*C) + r1*(1+r2/R)+r2;
a0 = r1/(R*C) + 1/C;
dx1dt = x(2);
dx2dt = (-a1 / a2 )* x(2) + (-a0/a2)*x(1) + b*Vs/a2 - (1+r2/R);% * dirac(t-0.01);
dxdt = [dx1dt ;dx2dt];
end
.

Community Treasure Hunt

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

Start Hunting!