0
$\begingroup$

I have solved a second order differential equation $\ddot R=f(t,R,\dot R)$ with a piecewise function $$P(t)=\begin{cases} P_0 & \text{ if } 0<t\leq40\times10^{-6} \\ \frac{(P_0+P_{min})}{2}+\frac{(P_0-P_{min})}{2}\cos(2\pi ft) & \text{ if } 40\times10^{-6}<t\leq60\times10^{-6} \\ P_{min} & \text{ if } 60\times10^{-6}<t\leq100\times10^{-6} \\ \frac{(P_0+P_{min})}{2}+\frac{(P_0-P_{min})}{2}\cos(2\pi ft) & \text{ if } 100\times10^{-6}<t\leq120\times10^{-6} \\ P_0 & \text{ if } t>120\times10^{-6} \end{cases}$$ enter image description here

I would like to plot, on the same graph with two $y$ axes (I guess with the plotyy command) both $t$ versus $R$ and $t$ versus $P(t)$. At the moment I am only able to plot $t$ versus $R$ plot(t,R(:,1)).

function overshoot()
R0 = 15e-6;
tspan = [0 160e-6];
[t,R] = ode45(@(t,R) DE(t,R,R0), tspan, [R0,0]);
t = t * 1e6;
R = R * 1e6;
plot(t,R(:,1))
end
%
function Rdot = DE(t,R,R0)
S = 0.073;
rho = 998;
mi = 1.005e-3; % kg / (m * s)
P0 = 101325;
Pvap = 2329.6;
Pmin = 1800;
f = 1e6 / 40;
%
Rdot = zeros(2,1);
Rdot(1) = R(2);
Rdot(2) = -1.5 * R(2) * R(2) / R(1) + 1 / (R(1) * rho) *...
    (Pvap - P(t,f,P0,Pmin) + (P0 - Pvap + 2 * S / R0) *...
    (R0 / R(1))^3 - 2 * S / R(1) - 4 * mi * R(2) / R(1));
end
%
function fval = P(t,f,P0,Pmin)
if (t <= 40e-6)
    fval = 101325;
elseif (t > 40e-6) && (t <= 60e-6)
    fval = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t);
elseif (t > 60e-6) && (t <= 100e-6)
    fval = 1800;
elseif (t > 100e-6) && (t <= 120e-6)
    fval = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t);
else
    fval = 101325;
end
end
$\endgroup$
  • $\begingroup$ Looks like this has been asked and answered on stackoverflow $\endgroup$ – Eric Kightley Dec 3 '16 at 15:57
  • $\begingroup$ I know that I have to use the plotyy command (which takes 4 arguments) but in this case I do not know how. The problem is with fval. For example: in this case, I am not even able to plot in a separeted graph t versum P. $\endgroup$ – Marco Dec 3 '16 at 15:59
  • 1
    $\begingroup$ This is a question purely about the correct syntax for a matlab script. We usually consider such questions off-topic here, and on topic for the matlab forums. $\endgroup$ – Wolfgang Bangerth Dec 4 '16 at 0:41
1
$\begingroup$

This can achieved by making the parameters P0 Pmin f global.

function overshoot()
global P0 Pmin f
R0 = 15e-6;

tspan = [0 160e-6]
[t,R] = ode45(@(t,R) DE(t,R,R0), tspan, [R0,0]);
%t = t * 1e6;
%R = R * 1e6;
figure(1)
plot(t,R(:,1),'r--')
hold on
plot(t,P(t,f,P0,Pmin),'o');
end  

function Rdot = DE(t,R,R0)
global P0 Pmin f
S = 0.073;
rho = 998;
mi = 1.005e-3; % kg / (m * s)
P0 = 101325;
Pvap = 2329.6;
Pmin = 1800;
f = 1e6 / 40;
%
Rdot = zeros(2,1);
Rdot(1) = R(2);
Rdot(2) = -1.5 * R(2) * R(2) / R(1) + 1 / (R(1) * rho) *...
    (Pvap - P(t,f,P0,Pmin) + (P0 - Pvap + 2 * S / R0) *...
    (R0 / R(1))^3 - 2 * S / R(1) - 4 * mi * R(2) / R(1));
end

function fval = P(t,f,P0,Pmin)
for i=1:length(t)
if (t(i) <= 40e-6)
    fval(i) = 101325;
elseif (t(i) > 40e-6) && (t(i) <= 60e-6)
    fval(i) = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t(i));
elseif (t(i) > 60e-6) && (t(i) <= 100e-6)
    fval(i) = 1800;
elseif (t(i) > 100e-6) && (t(i) <= 120e-6)
    fval(i) = (P0 + Pmin) / 2 + (P0 - Pmin) / 2 * cos(2 * pi * f * t(i));
else
    fval(i) = 101325;
end
end    
end

The output plot which now has two $y$ axes, i.e, $R(t)$ and $P(t)$. But the plot doesn't make any sense because the maximum for $R(t)$ is approximately $8\times10^{-5}$ and for $P(t)$ is $\approx 10\times10^{4}$.

Anyways the desired combined plot looks like this enter image description here

I hope this helps.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.