计算极限环内的不动点

计算科学 固定点
2021-12-02 23:50:26

所以我正在使用一个相当复杂的动态系统。而不是全部写出来。如果您只是克隆我的 git 存储库,可能会更容易。

git clone https://gitlab.com/mdornfe1/vortex_acoustic

系统的动力学被编码在 vortex_acoustic.py 中的函数流中。如果我从休息中整合系统,我会发现某些参数值的振荡。该系统是 13 维的,但这里是第一个变量与时间的关系图 在此处输入图像描述

所以看起来系统正在接近极限环吸引子。根据我对动力系统的理解,该极限环内应该有一个不稳定的不动点。我想尝试计算它,以便分析它的稳定性。我的想法是固定点应该接近振荡量的平均值,因此我将其用作数值根求解器的起始值。这是查找固定点的代码,但您也可以只运行 stack_exchange.py。

import pickle, numpy as np, matplotlib.pyplot as plt
import vortex_acoustic as va
from scipy import optimize
from constants import *

T = 5000
dt = 0.05
tn_transients = int(0.9 * T / dt)

def calc_fixed_point(seed, params):
    return optimize.fsolve(va.flow_star, seed, params)

def calc_closest_fixed_point(y, params):
    seed = np.mean(y[tn_transients:, :], 0)
    seed[Nq+1:2*Nq+1] = 0
    y_star = calc_fixed_point(seed, params)

    return y_star, seed

y = np.load('simulation.npy')
with open('params.p', 'r') as f:
    params = pickle.load(f)

plt.plot(y[:,0])

y_star, seed = calc_closest_fixed_point(y, params)

函数 calc_closest_fixed_point 获取模拟的输出和传递给模拟的参数。它计算瞬态消失后波动量的平均值。它将导数设置为 0,并使用该数量种子作为 fsolve 的起点。然后函数 calc_fixed_point 找到系统 y_star 的不动点。我认为这种方法会表明种子和 y_star 彼此非常接近。但是对于这个例子,我看到它们完全不同。此外,实际值不在极限循环内。

print(seed)
>>>array([ 0.18307736, -0.11207637,  0.00382286,  0.        ,  0.        ,
       -0.02016129,  0.08881281, -0.07852504,  0.07439292,  0.0109214 ,
       -0.00894938, -0.02896082,  0.00833759])
print(y_star)
>>>array([  1.06466138e+00,  -6.65605305e-01,   2.15920825e-02,
        -5.18569988e-25,   1.44087666e-24,   7.55288134e-03,
         5.17875853e-01,  -4.41152924e-01,   4.72080951e-01,
         4.86813549e-02,  -6.04678149e-02,  -1.42533599e-01,
         5.90028044e-02])

种子也很明显不是一个固定点。如果您在种子时计算函数流,则导数与 0 相差甚远。

va.flow(seed,0,*params)
>>>array([-0.00543102,  0.        ,  0.        , -0.00113827, -0.00042674,
       -0.00973662, -0.00775962,  0.01449168,  0.01369892,  0.00331889,
        0.00502761, -0.00026469, -0.00076228])

有谁知道发生了什么?我假设在极限循环内有一个固定点是正确的,对吗?

0个回答
没有发现任何回复~