求解和绘制分段 ODE答案

【问题标题】：Solving and plotting a piecewise ODE求解和绘制分段 ODE
【发布时间】：2020-12-20 07:14:53
【问题描述】：

我有一个函数dφ/dt = γ - F(φ)（其中F(φ)——a 是2π-周期函数）和函数图F(φ)。

我需要创建一个程序，为 γ 的不同值（γ = 0.1、0.5、0.95、1.05、2、5）输出 6 个 φ(t) 的图，和t∈[0,100]。

这里是F(φ)函数的定义：

      -φ/a - π/a,    if φ ∈ [-π, -π + a]
      -1,            if φ ∈ [-π + a, - a] 
F(φ) = φ/a,          if φ ∈ [- a, a]
       1,            if φ ∈ [a, π - a]
      -φ/a + π/a,    if φ ∈ [π - a, π]

                 ^ F(φ)
                 |
                 |1   ______
                 |   /|     \
                 |  / |      \
                 | /  |       \      φ
__-π_______-a____|/___|________\π____>
   \        |   /|0    a
    \       |  / |
     \      | /  |
      \     |/   |
       ¯¯¯¯¯¯    |-1

我的问题是我不知道在边界和初始条件方面给ode45 提供什么输入。我所知道的是φ(t) 的演变必须是连续的。

这是γ = 0.1案例的代码：

hold on;
df1dt = @(t,f1) 0.1 - f1 - 3.14;
df2dt = @(t,f2)- 1;
df3dt = @(t,f3) 0.1 + f3;
df4dt = @(t,f4)+1;
df5dt = @(t,f5) 0.1 - f5 + 3.14;
[T1,Y1] = ode45(df1dt, ...);
[T2,Y2] = ode45(df2dt, ...);
[T3,Y3] = ode45(df3dt, ...);
[T4,Y4] = ode45(df4dt, ...);
[T5,Y5] = ode45(df5dt, ...);
plot(T1,Y1);
plot(T2,Y2);
plot(T3,Y3);
plot(T4,Y4);
plot(T5,Y5);
hold off; 
title('\gamma = 0.1')

【问题讨论】：

通常，在求解一阶 ODE 时，需要 1 个初始条件，如果没有这些初始条件，则问题是欠定义的。在这种情况下，您可以猜测一个初始条件，但如果猜测错误，解决方案将偏离一个常数。
与math.stackexchange.com/questions/3955811/…交叉发布

标签： matlab plot ode piecewise periodicity

【解决方案1】：

我们先定义F(φ,a)：

function out = F(p, a)
phi = mod(p,2*pi);
out = (0      <= phi & phi < a     ).*(phi/a) ...
    + (a      <= phi & phi < pi-a  ).*(1) ...
    + (pi-a   <= phi & phi < pi+a  ).*(-phi/a + pi/a) ...
    + (pi+a   <= phi & phi < 2*pi-a).*(-1) ...
    + (2*pi-a <= phi & phi < 2*pi  ).*(phi/a - 2*pi/a);
end

对于某些示例输入给出：

使用绘图代码：

x = linspace(-3*pi, 3*pi, 200);
a = pi/6;

figure(); plot(x,F(x, a));
xlim([-3*pi,3*pi]);
xticks(-3*pi:pi:3*pi);
xticklabels((-3:3)+ "\pi");
grid on; grid minor
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):pi/6:ax.XAxis.Limits(2);

从那里您无需再担心范围，只需调用ode45：

% Preparations:
a = pi/6;
g = [0.1, 0.5, 0.95, 1.05, 2, 5]; % γ
phi0 = 0; % you need to specify the correct initial condition (!)
tStart = 0;
tEnd = 100;
% Calling the solver:
[t, phi] = arrayfun(@(x)ode45(@(t,p)x-F(p,a), [tStart, tEnd], phi0), g, 'UniformOutput', false);
% Plotting:
plotData = [t; phi];
figure(); plot(plotData{:});
legend("γ=" + g, 'Location', 'northwest');

导致：

【讨论】：