求解并绘制分段 ODE

Question

我有一个函数d\xcf\x86/dt = \xce\xb3 - F(\xcf\x86)（其中F(\xcf\x86)-- a 是2\xcf\x80周期函数）和该函数的图形F(\xcf\x86)。

\n

我需要创建一个程序，输出 6 个\xcf\x86(t)不同值\xce\xb3( \xce\xb3 = 0.1、0.5、0.95、1.05、2、5)和 的图t\xe2\x88\x88[0,100]。

\n

这是函数的定义F(\xcf\x86)：

\n

      -\xcf\x86/a - \xcf\x80/a,    if \xcf\x86 \xe2\x88\x88 [-\xcf\x80, -\xcf\x80 + a]\n      -1,            if \xcf\x86 \xe2\x88\x88 [-\xcf\x80 + a, - a] \nF(\xcf\x86) = \xcf\x86/a,          if \xcf\x86 \xe2\x88\x88 [- a, a]\n       1,            if \xcf\x86 \xe2\x88\x88 [a, \xcf\x80 - a]\n      -\xcf\x86/a + \xcf\x80/a,    if \xcf\x86 \xe2\x88\x88 [\xcf\x80 - a, \xcf\x80]\n

\n

                 ^ F(\xcf\x86)\n                 |\n                 |1   ______\n                 |   /|     \\\n                 |  / |      \\\n                 | /  |       \\      \xcf\x86\n__-\xcf\x80_______-a____|/___|________\\\xcf\x80____>\n   \\        |   /|0    a\n    \\       |  / |\n     \\      | /  |\n      \\     |/   |\n       \xc2\xaf\xc2\xaf\xc2\xaf\xc2\xaf\xc2\xaf\xc2\xaf    |-1\n

\n

我的问题是我不知道ode45在边界和初始条件方面要给出什么输入。我所知道的是，进化 \xcf\x86(t)必须是连续的。

\n

这是以下情况的代码\xce\xb3 = 0.1：

\n

hold on;\ndf1dt = @(t,f1) 0.1 - f1 - 3.14;\ndf2dt = @(t,f2)- 1;\ndf3dt = @(t,f3) 0.1 + f3;\ndf4dt = @(t,f4)+1;\ndf5dt = @(t,f5) 0.1 - f5 + 3.14;\n[T1,Y1] = ode45(df1dt, ...);\n[T2,Y2] = ode45(df2dt, ...);\n[T3,Y3] = ode45(df3dt, ...);\n[T4,Y4] = ode45(df4dt, ...);\n[T5,Y5] = ode45(df5dt, ...);\nplot(T1,Y1);\nplot(T2,Y2);\nplot(T3,Y3);\nplot(T4,Y4);\nplot(T5,Y5);\nhold off; \ntitle(\'\\gamma = 0.1\')\n

\n

Answer 1

让我们首先定义F(\xcf\x86,a)：

\n

function out = F(p, a)\nphi = mod(p,2*pi);\nout = (0      <= phi & phi < a     ).*(phi/a) ...\n    + (a      <= phi & phi < pi-a  ).*(1) ...\n    + (pi-a   <= phi & phi < pi+a  ).*(-phi/a + pi/a) ...\n    + (pi+a   <= phi & phi < 2*pi-a).*(-1) ...\n    + (2*pi-a <= phi & phi < 2*pi  ).*(phi/a - 2*pi/a);\nend\n

\n

对于某些示例输入，给出：\n

\n

使用绘图代码：

\n

x = linspace(-3*pi, 3*pi, 200);\na = pi/6;\n\nfigure(); plot(x,F(x, a));\nxlim([-3*pi,3*pi]);\nxticks(-3*pi:pi:3*pi);\nxticklabels((-3:3)+ "\\pi");\ngrid on; grid minor\nax = gca;\nax.XAxis.MinorTick = \'on\';\nax.XAxis.MinorTickValues = ax.XAxis.Limits(1):pi/6:ax.XAxis.Limits(2);\n

\n

从那里你不需要再担心范围，只需调用ode45：

\n

% Preparations:\na = pi/6;\ng = [0.1, 0.5, 0.95, 1.05, 2, 5]; % \xce\xb3\nphi0 = 0; % you need to specify the correct initial condition (!)\ntStart = 0;\ntEnd = 100;\n% Calling the solver:\n[t, phi] = arrayfun(@(x)ode45(@(t,p)x-F(p,a), [tStart, tEnd], phi0), g, \'UniformOutput\', false);\n% Plotting:\nplotData = [t; phi];\nfigure(); plot(plotData{:});\nlegend("\xce\xb3=" + g, \'Location\', \'northwest\');\n

\n

导致：

\n

\n