How to solve Piece wise Linear ODE using ODE45

Question

0 voti

Dear Friends, I hope you guys are doing good. I have small doubt on "How to solve Piecewise Linear ODE using ODE45. I am hereby attaching my code (Basically Energy Harvesting). Please kindly look at my code. I want to plot Frequency vs Volts.

          function xdot = Linear2DOF(t,x)
          global m1 m2 k1 k2 c1 c2 W theta R cs K d p1 p2
          if((x(1)-x(3))<-d)
              p1 = (K/m1)*(x(1)-x(3)+d);
              p2 = (K/m2)*(x(1)-x(3)+d);
          elseif((x(1)-x(3))>d)
              p1 = (K/m1)*(x(1)-x(3)-d);
              p2 = (K/m2)*(x(1)-x(3)-d);
          else
              p1 = 0;
              p2 = 0;
          end
              xdot = [x(2); -((c1/m1)*x(2))-((k1/m1)*x(1))-((c2/m1)*(x(2)-x(4)))-((k2/m1)*(x(1)-x(3)))-((theta/m1)*x(5))-p1+(2*cos(W*t)); x(4); -((c2/m2)*(x(4)-x(2)))-((k2/m2)*(x(3)-x(1)))+p2+(2*cos(W*t)); -(x(5)/(R*cs))+((theta*x(2))/cs)]; 
          end

My script file is given below

clc;
clear all;
close all;
global m1 m2 k1 k2 c1 c2 W theta R cs K d 
m1    = 0.056;
m2    = 0.0084;
k1    = 1500;
k2    = 144;
c1    = 0.1833;
c2    = 0.0132;
theta = 1.3*10^-3;
R     = 10^20;
cs    = 180*10^-9;
K    = 15000;
d     = 0.005;
omega = [];
Disp1 = [];
Velo1 = [];
Disp2 = [];
Velo2 = [];
Volts = [];
x0    = [0; 0; 0; 0; 0];
tspan = [0 20];
for W = 100:10:300
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,x] = ode23(@Linear2DOF,tspan,x0);
omega(end+1) = W/6.28; % Frequency in Hz
Disp1(end+1) = max(x(:,1));
Velo1(end+1) = max(x(:,2));
Disp2(end+1) = max(x(:,3));
Velo2(end+1) = max(x(:,4));
Volts(end+1) = max(x(:,5));
end