Чувствительность уравнения Дуффинга 

Изучение чувствительности решений уравнения Дуффинга к изменениям в начальных условиях.

Для упрощения последующих вычислений используем команду ParametricNDSolveValue, в которой начальные условия являются параметрами, а остальные параметры фиксированы.

In[1]:=

X xp = ParametricNDSolveValue[{x''[ t] + \[Delta] x'[t] + \[Alpha] x[ t] + \[Beta] x[t]^3 == \[Gamma] Cos[\[Omega] t], x[0] == a, x'[0] == b} /. {\[Delta] -> 0.15, \[Alpha] -> -1, \[Beta] -> 1, \[Gamma] -> 0.3, \[Omega] -> 1}, x, {t, 0, 100}, {a, b}]

Out[1]=

Производные объекта ParametricFunction по параметрам a и b являются чувствительностями относительно начальных значений для x и его производной x соответственно.

In[2]:=

X {xpa = Head[D[xp[a, b], a]], xpb = Head[D[xp[a, b], b]]}

Out[2]=

Отображение графика решения, с начальными значениями {0, 0} и графиков чувствительностей в логарифмическом масштабе.

In[3]:=

X GraphicsRow[{Plot[xp[0, 0][t], {t, 0, 100}], LogPlot[Abs[xpa[0, 0][t]], {t, 0, 100}], LogPlot[Abs[xpb[0, 0][t]], {t, 0, 100}]}]

Out[3]=

Огромные значения чувствительностей указывают на то, что даже малейшее изменение в одном из параметров повлечёт существенное отклонение в решении. Ниже приведены графики решений с малым отклонением от нулевых начальных условий по каждому из направлений.

In[4]:=

X \[Epsilon] = 10^-6; Plot[{xp[0, 0][t], xp[\[Epsilon], 0][t], xp[0, \[Epsilon]][t]}, {t, 0, 100}, ImageSize -> Medium]

Out[4]=

Ещё один способ визуализации чувствительности - это произвести её локальное разбиение на составляющие, параллельные и перпендикулярные к направлению траектории на фазовой плоскости.

Определим функцию, находящую вектор, перпендикулярный к траектории как функцию от времени.

In[5]:=

X perp[a_, b_][t_?NumberQ] := Module[{f = xp[a, b], dt, do}, dt = Normalize[{f'[t], f''[t]}]; do = Reverse[dt]*{1, -1} ];

Теперь определим функции, дающие модуль компонент чувствительности по отношению к начальному значению x в параллельном и перпендикулярном направлениях.

In[6]:=

X spar[a_, b_][t_?NumberQ] := Module[{s = xpa[a, b], f = xp[a, b]}, Abs[Normalize[{f'[t], f''[t]}].{s[t], s'[t]}]]

In[7]:=

X sperp[a_, b_][t_?NumberQ] := Module[{s = xpa[a, b]}, Abs[perp[a, b][t].{s[t], s'[t]}]]

Для этого уравнения эти две компоненты растут со сравнимыми скоростями.

Разбиение на компоненты позволяет строить занятные визуализации на фазовой плоскости. Изменение масштаба чувствительности позволяет наблюдать её эффект для различных временных интервалов.

In[8]:=

X scale = 10^-8; Show[ ParametricPlot[{xp[0, 0][t], xp[0, 0]'[t]} + s sperp[0, 0][t] perp[0, 0][t], {t, 0, 100}, {s, -scale, scale}, Mesh -> None, PlotPoints -> {1000, 2}], ParametricPlot[{xp[0, 0][t], xp[0, 0]'[t]}, {t, 0, 25}, PlotStyle -> Red], ImageSize -> Medium]

Out[8]=