求解和绘制分段 ODE

问题描述

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

我需要创建一个程序，针对不同的 φ(t) (γ,γ = 0.1,0.5,{{1} }}、0.95、1.05) 和 2。

这里是 5 函数的定义：

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)

解决方法

让我们先定义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');

结果：