计算科学 - Matlab 中的 ODE：如何绘制结果和分段函数 - 吾爱随笔录

Matlab 中的 ODE：如何绘制结果和分段函数

计算科学 matlab

2021-12-11 17:57:23

我已经用分段函数求解了二阶微分方程 $\ddot R=f(t,R,\dot R)$

P (t) = {\begin{cases} P_{0} & if 0 < t \leq 40 \times 10^{- 6} \\ \frac{(P_{0} + P_{m i n})}{2} + \frac{(P_{0} - P_{m i n})}{2} \cos (2 π f t) & if 40 \times 10^{- 6} < t \leq 60 \times 10^{- 6} \\ P_{m i n} & if 60 \times 10^{- 6} < t \leq 100 \times 10^{- 6} \\ \frac{(P_{0} + P_{m i n})}{2} + \frac{(P_{0} - P_{m i n})}{2} \cos (2 π f t) & if 100 \times 10^{- 6} < t \leq 120 \times 10^{- 6} \\ P_{0} & if t > 120 \times 10^{- 6} \end{cases}

$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}$

我想用两个轴（我猜是命令）与和与。目前我只能绘制与的图。 $y$ plotyy $t$ $R$ $t$ $P(t)$ $t$ $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

1个回答

这可以通过使参数P0 Pmin f全局化来实现。

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

现在有两个轴的输出图，即和。但是该图没有任何意义，因为的最大值约为而的。 $y$ $R(t)$ $P(t)$ $R(t)$ $8\times10^{-5}$ $P(t)$ $\approx 10\times10^{4}$

无论如何，所需的组合情节看起来像这样

我希望这有帮助。

其它你可能感兴趣的问题

上一篇两级迭代器 C++ 下一篇非厄米系数矩阵在自伴算子离散化中的作用