Как я могу заставить свое решение сходиться для определения фиксированной точки для нелинейного набора уравнений? [закры

Как я могу заставить свое решение сходиться для определения фиксированной точки для нелинейного набора уравнений? [закры ⇐ Python

1 сообщение • Страница 1 из 1

Anonymous

Как я могу заставить свое решение сходиться для определения фиксированной точки для нелинейного набора уравнений? [закры

Цитата

Сообщение Anonymous » 20 ноя 2025, 08:48

В связи с моим предыдущим вопросом, в котором я упомянул, что пытаюсь написать код на Python, используя Sympy, с целью определить фиксированную точку нелинейной системы. Я импортирую константы из файла .mat.

Код не может определить фиксированные точки и не сходится. Сначала я пытался использовать для поиска Sympy.solve, но это не сработало. Затем я попробовал numpy.nsolve, который тоже не сработал. Затем я выбрал корневую функцию, которая тоже не сработала.
Функция Sympy sp.Abs() для поиска абсолютного значения также вообще не работает при поиске того же самого, хотя g1 является символической.

import Sympy как sp
из Sympy import nsolve, Abs, символы
import numpy как np
из scipy.optimize import root
из pymatreader import read_mat
импортировать matplotlib.pyplot как plt
импортировать pandas как pd
импортировать трассировку
импортировать математические данные
импортировать предупреждения
warnings.filterwarnings("ignore", Category=RuntimeWarning)

--- Загрузить файл MAT ---

file_path = r'E:\RC\data&plot\data1.mat'

data = read_mat(file_path)
data1 = data['data1']
vars_ = data1['variable']
consts = data1['constant']

results = [] # сохраняем решения для каждого столбца
i=9998

num_cols = vars_['Q'].shape1
g1_fp1 = [Нет] * num_cols
g1_fp2 = [Нет] * num_cols
g2_fp = [Нет] * num_cols
x_val=[]
J=[]
fp=[]
g1_mean = []
g1_diff = []
pad_list=[]
rows=[]
eig_list=[]
eig_vals_real=[]

Символические параметры
x, g1, g2, Q, Gam, I_expr, den, chi0, chi1, D0, D1, Alpha_chi, альфа_D, chi_, beta_, lam, mu_c, > M2, M3, M4, M5 = sp.symbols(
'x, g1, g2, Q, Gam, I, den, chi0, chi1, D0, D1, альфа_chi, Alpha_D, chi_, beta_, lam, mu_c, > M2, M3, M4, M5', real=True
)

Построение символьных вспомогательных выражений
A = (-1/den2) * (Gam * chi1) / ( (Q * D1) + g1 * (D0 * chi1 - D1 * chi0))
B = (-1/den2) * (Q * D1) / ( (Gam * chi1) + g2 * (chi0 * D1 - D0 * chi1) )

v_eQ = (-1 / den2) * (Gam * chi1) / ( (Q * D1 + g1 * (D0 * chi1 - D1 * chi0)) * g12 )
v_eGam = (-1 / den2) * (Q * D1) / ( (Gam * chi1 + g2 * (chi0 * D1 - D0 * chi1)) * g22 )

Уравнения (символические)
f1 = Q - (chi0 + chi1 * I_expr / (1 + Alpha_chi * v_eQ2)) * g1
f2 = Gam - (D0 + D1 * I_expr / (1 + Alpha_D * v_eGam2)) * g2

abs_term = sp.Abs(-g1 - mu_c_val)
f3 = (chi_ / beta_) * (abs_term - lam * A2 * g14)-I_expr

def root_to_float(r, tol_imag=1e-9):
попробуйте:
rv = complex(sp.N(r))
if abs(rv.imag) < tol_imag:
return float(rv.real)
кроме исключения:
pass
return np.nan

for j in range(num_cols):
print(f"\n--- Столбец {j} ---")

--- Извлечь значения скалярного столбца -----
def extract_scalar(val, j):
"""Вернуть val[j], если val похоже на массив, в противном случае вернуть val как число с плавающей запятой."""
if isinstance(val, np.ndarray):
if val.ndim == 0:
return float(val) # скаляр ndarray
elif val.ndim == 1:
return float(val[j]) # 1D массив
elif val.ndim == 2:
return float(val[:, j]) # столбец из 2D массива
else:
raise ValueError(f"Неподдерживаемая форма для переменной: {val.shape}")
else:
return float(val)

x_val = extract_scalar(np.array(vars_['x']), j)
Q_val = float(vars_['Q'][0, j])
Gam_val = float(vars_['Gamma'][0, j])
I_val = float(vars_['turbulence'][0, j])
den_val = float(vars_['density'][0, j])
g1_val = float(vars_['gradpressure'][0, j])
g2_val = float(vars_['graddensity'][0, j])
x_val = float(vars_['x'][j]) # x — это 1D, поэтому только j-й элемент

--- Константы ---
consts = data['data1']['constant']

chi0_val = consts.get('chi0')
chi1_val = consts.get('chi1')
D0_val = consts.get('D0')
D1_val = consts.get('D1')
alpha_chi_val = consts.get('alphachi')
alpha_D_val = consts.get('alphaD')
chi_val = consts.get('chi_growth', consts.get('chi_', None)) # попробуйте общие имена
beta_val = consts.get('beta_suppression', consts.get('beta_', None))
lam_val = consts.get('lambda_suppress', consts.get('lambda', None))
mu_c_val = float(consts.get('critgradpressure',None))

--- Производное ---
M2_val = (Gam_val2) * (chi1_val2)
M3_val = (Q_val2) * (D1_val2)
M4_val = (D0_valchi1_val - chi0_valD1_val)2
M5_val = 2(D0_valchi1_val - chi0_val*D1_val)

--- Подстановка ---
subs_ = {
x: x_val, Q: Q_val, Gam: Gam_val, I_expr: I_val, den: den_val,
chi0: chi0_val, chi1: chi1_val, D0: D0_val, D1: D1_val,
alpha_chi: Alpha_chi_val, Alpha_D: Alpha_D_val,
/>chi_: chi_val, beta_: beta_val,
lam: lam_val, mu_c: mu_c_val,
M2: M2_val, M3: M3_val, M4: M4_val, M5: M5_val

def f1_numeric(g1, Q_val, Gam_val, M2, M3, M4, M5):
# Ваш числовой эквивалент f1
return (Q - (chi0 + chi1 * I_expr / (1 + Alpha_chi * v_eQ**2)) * g1)

def f2_numeric(g2, Q_val, Gam_val, M2, M3, M4, M5):
возврат (Gam - (D0 + D1 * I_expr / (1 + альфа_D * v_eGam2)) * g2)
def f3_numeric(g1, chi_val, beta_val, lam_val, A, mu_c_val, I_val):
возврат (chi_val / beta_val) *

(np.abs(-g1 - mu_c_val) - lam_val * A2 * g1**4) - I_val

--- Найти фиксированные точки ---
g1_guess = 0,031
g2_guess = 2,3

попробуйте:
#r3 = sp.solve(f3_num, g1)
r3 = root(f3_num, g1, g1_guess)
if r3:
g1_fp1[j] = root_to_float(r3[0])
кроме исключения:
g1_fp1[j] = np.nan
print("g1_fp1=",g1_fp1)

try:
#r1 = sp.solve(f1_num, g1)
r1 = root(f1_num, g1, g1_guess)
if r1:
g1_fp2[j] = root_to_float(r1[0])
кроме исключения:
g1_fp2[j] = np.nan
print("g1_fp2=",g1_fp2)

попробуй:
#r2 = sp.solve(f2_num, g2)
r2 = root(f2_num, g2, g2_guess)
if r2:
g2_fp[j] = root_to_float(r2[0])
кроме исключения:
g2_fp[j] = np.nan
print("g2_fp=",g2_fp)

g1_fp1_arr = np.array(g1_fp1, dtype=float)
g1_fp2_arr = np.array(g1_fp2, dtype=float)
g2_fp_arr = np.array(g2_fp, dtype=float)

если g1_fp1_arr не имеет значения None и g1_fp2_arr не имеет значения None:
mean_j = np.nanmean([g1_fp1_arr, g1_fp2_arr])
diff_j = g1_fp1_arr - g1_fp2_arr
else:
mean_j = np.nan
diff_j = np.nan
vals = [g1_fp1_arr, g1_fp2_arr]

g1_mean = np.nanmean(np.vstack([g1_fp1_arr, g1_fp2_arr]), axis=0)
g1_diff = g1_fp1_arr - g1_fp2_arr

defpad_list(lst, length):
return lst + [np.nan] * (длина - len(lst))

#Преобразовать массив Numpy в список перед заполнением
g1_fp1_arr = Pad_list(g1_fp1_arr.tolist(), max_len)
g1_fp2_arr = Pad_list(g1_fp2_arr.tolist(), max_len)
g1_mean = Pad_list(g1_mean.tolist(), max_len)
g1_diff = Pad_list(g1_diff.tolist(), max_len)
g2_fp_arr = Pad_list(g2_fp_arr.tolist(), max_len)

output_path = r"E:\RC\fp_values.xlsx"

x_arr = np.array(vars_['x']). Flatten()

Определите якобиан символически
J = sp.Matrix([
[sp.diff(f1, g1), sp.diff(f1, g2), sp.diff(f1, I_expr)],
[sp.diff(f2, g1), sp.diff(f2, g2), sp.diff(f2, I_expr)],
[sp.diff(f3, g1), sp.diff(f3, g2), sp.diff(f3, I_expr)]
])
sp.pprint(J)

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

Подробнее здесь: https://stackoverflow.com/questions/798 ... -nonlinear

1763617702

Anonymous

В связи с моим предыдущим вопросом, в котором я упомянул, что пытаюсь написать код на Python, используя Sympy, с целью определить фиксированную точку нелинейной системы. Я импортирую константы из файла .mat.
[list]
[*]Код не может определить фиксированные точки и не сходится. Сначала я пытался использовать для поиска Sympy.solve, но это не сработало. Затем я попробовал numpy.nsolve, который тоже не сработал.  Затем я выбрал корневую функцию, которая тоже не сработала.

[*]Функция Sympy sp.Abs() для поиска абсолютного значения также вообще не работает при поиске того же самого, хотя g1 является символической.

[/list]

import Sympy как sp
из Sympy import nsolve, Abs, символы
import numpy как np
из scipy.optimize import root
из pymatreader import read_mat
импортировать matplotlib.pyplot как plt
импортировать pandas как pd
импортировать трассировку
импортировать математические данные
импортировать предупреждения
warnings.filterwarnings("ignore", Category=RuntimeWarning)

--- Загрузить файл MAT ---

file_path = r'E:\RC\data&plot\data1.mat'


data = read_mat(file_path)
data1 = data['data1']
vars_ = data1['variable']
consts = data1['constant']


results = [] # сохраняем решения для каждого столбца
i=9998

num_cols = vars_['Q'].shape1
g1_fp1 = [Нет] * num_cols
g1_fp2 = [Нет] * num_cols
g2_fp = [Нет] * num_cols
x_val=[]
J=[]
fp=[]
g1_mean = []
g1_diff = []
pad_list=[]
rows=[]
eig_list=[]
eig_vals_real=[]


Символические параметры
x, g1, g2, Q, Gam, I_expr, den, chi0, chi1, D0, D1, Alpha_chi, альфа_D, chi_, beta_, lam, mu_c, > M2, M3, M4, M5 = sp.symbols(
'x, g1, g2, Q, Gam, I, den, chi0, chi1, D0, D1, альфа_chi, Alpha_D, chi_, beta_, lam, mu_c, > M2, M3, M4, M5', real=True
)


Построение символьных вспомогательных выражений
A = (-1/den[b]2) * (Gam * chi1) / ( (Q * D1) + g1 * (D0 * chi1 - D1 * chi0))
B = (-1/den[/b]2) * (Q * D1) / ( (Gam * chi1) + g2 * (chi0 * D1 - D0 * chi1) )


v_eQ = (-1 / den[b]2) * (Gam * chi1) / ( (Q * D1 + g1 * (D0 * chi1 - D1 * chi0)) * g1[/b]2 )
v_eGam = (-1 / den[b]2) * (Q * D1) / ( (Gam * chi1 + g2 * (chi0 * D1 - D0 * chi1)) * g2[/b]2 )


Уравнения (символические)
f1 = Q - (chi0 + chi1 * I_expr / (1 + Alpha_chi * v_eQ[b]2)) * g1
f2 = Gam - (D0 + D1 * I_expr / (1 + Alpha_D * v_eGam[/b]2)) * g2


abs_term = sp.Abs(-g1 - mu_c_val)
f3 = (chi_ / beta_) * (abs_term - lam * A[b]2 * g1[/b]4)-I_expr


def root_to_float(r, tol_imag=1e-9):
попробуйте:
rv = complex(sp.N(r))
if abs(rv.imag) < tol_imag:
return float(rv.real)
кроме исключения:
pass
return np.nan


for j in range(num_cols):
print(f"\n--- Столбец {j} ---")


--- Извлечь значения скалярного столбца -----
def extract_scalar(val, j):
"""Вернуть val[j], если val похоже на массив, в противном случае вернуть val как число с плавающей запятой."""
if isinstance(val, np.ndarray):
if val.ndim == 0:
return float(val) # скаляр ndarray
elif val.ndim == 1:
return float(val[j]) # 1D массив
elif val.ndim == 2:
return float(val[:, j]) # столбец из 2D массива
else:
raise ValueError(f"Неподдерживаемая форма для переменной:  {val.shape}")
else:
return float(val)


x_val = extract_scalar(np.array(vars_['x']), j)
Q_val = float(vars_['Q'][0, j])
Gam_val = float(vars_['Gamma'][0, j])
I_val = float(vars_['turbulence'][0, j])
den_val = float(vars_['density'][0, j])
g1_val = float(vars_['gradpressure'][0, j])
g2_val = float(vars_['graddensity'][0, j])
x_val = float(vars_['x'][j]) # x — это 1D, поэтому только j-й элемент


--- Константы ---
consts = data['data1']['constant']


chi0_val = consts.get('chi0')
chi1_val = consts.get('chi1')
D0_val = consts.get('D0')
D1_val = consts.get('D1')
alpha_chi_val = consts.get('alphachi')
alpha_D_val = consts.get('alphaD')
chi_val = consts.get('chi_growth', consts.get('chi_', None)) # попробуйте общие имена
beta_val = consts.get('beta_suppression', consts.get('beta_', None))
lam_val = consts.get('lambda_suppress', consts.get('lambda', None))
mu_c_val = float(consts.get('critgradpressure',None))


--- Производное ---
M2_val = (Gam_val[b]2) * (chi1_val[/b]2)
M3_val = (Q_val[b]2) * (D1_val[/b]2)
M4_val = (D0_valchi1_val - chi0_valD1_val)2
M5_val = 2(D0_valchi1_val - chi0_val*D1_val)


--- Подстановка ---
subs_ = {
x: x_val, Q: Q_val, Gam: Gam_val, I_expr: I_val, den: den_val,
chi0: chi0_val, chi1: chi1_val, D0: D0_val, D1: D1_val,
alpha_chi: Alpha_chi_val, Alpha_D: Alpha_D_val,
/>chi_: chi_val, beta_: beta_val,
lam: lam_val, mu_c: mu_c_val,
M2: M2_val, M3: M3_val, M4: M4_val, M5: M5_val



def f1_numeric(g1, Q_val, Gam_val, M2, M3, M4, M5):
# Ваш числовой эквивалент f1
return (Q - (chi0 + chi1 * I_expr / (1 + Alpha_chi * v_eQ**2)) * g1)


def f2_numeric(g2, Q_val, Gam_val, M2, M3, M4, M5):
возврат (Gam - (D0 + D1 * I_expr / (1 + альфа_D * v_eGam[b]2)) * g2)
def f3_numeric(g1, chi_val, beta_val, lam_val, A, mu_c_val, I_val):
возврат (chi_val / beta_val) *

(np.abs(-g1 - mu_c_val) - lam_val * A[/b]2 * g1**4) - I_val


--- Найти фиксированные точки ---
g1_guess = 0,031
g2_guess = 2,3


попробуйте:
#r3 = sp.solve(f3_num, g1)
r3 = root(f3_num, g1, g1_guess)
if r3:
g1_fp1[j] = root_to_float(r3[0])
кроме исключения:
g1_fp1[j] = np.nan
print("g1_fp1=",g1_fp1)


try:
#r1 = sp.solve(f1_num, g1)
r1 = root(f1_num, g1, g1_guess)
if r1:
g1_fp2[j] = root_to_float(r1[0])
кроме исключения:
g1_fp2[j] = np.nan
print("g1_fp2=",g1_fp2)


попробуй:
#r2 = sp.solve(f2_num, g2)
r2 = root(f2_num, g2, g2_guess)
if r2:
g2_fp[j] = root_to_float(r2[0])
кроме исключения:
g2_fp[j] = np.nan
print("g2_fp=",g2_fp)


g1_fp1_arr = np.array(g1_fp1, dtype=float)
g1_fp2_arr = np.array(g1_fp2, dtype=float)
g2_fp_arr = np.array(g2_fp, dtype=float)


если g1_fp1_arr не имеет значения None и g1_fp2_arr не имеет значения None:
mean_j = np.nanmean([g1_fp1_arr, g1_fp2_arr])
diff_j = g1_fp1_arr - g1_fp2_arr
else:
mean_j = np.nan
diff_j = np.nan
vals = [g1_fp1_arr, g1_fp2_arr]


g1_mean = np.nanmean(np.vstack([g1_fp1_arr, g1_fp2_arr]), axis=0)
g1_diff = g1_fp1_arr - g1_fp2_arr


defpad_list(lst, length):
return lst + [np.nan] * (длина - len(lst))


#Преобразовать массив Numpy в список перед заполнением
g1_fp1_arr = Pad_list(g1_fp1_arr.tolist(), max_len)
g1_fp2_arr = Pad_list(g1_fp2_arr.tolist(), max_len)
g1_mean = Pad_list(g1_mean.tolist(), max_len)
g1_diff = Pad_list(g1_diff.tolist(), max_len)
g2_fp_arr = Pad_list(g2_fp_arr.tolist(), max_len)


output_path = r"E:\RC\fp_values.xlsx"


x_arr = np.array(vars_['x']).  Flatten()


Определите якобиан символически
J = sp.Matrix([
[sp.diff(f1, g1), sp.diff(f1, g2), sp.diff(f1, I_expr)],
[sp.diff(f2, g1), sp.diff(f2, g2), sp.diff(f2, I_expr)],
[sp.diff(f3, g1), sp.diff(f3, g2), sp.diff(f3, I_expr)]
])
sp.pprint(J)

[img]https://i.sstatic.net/1Kw5Gg63.png[/img]
Выходные данные, которые я получаю, показаны на прикрепленном изображении как Нет.
ваши рекомендации по правильной оценке фиксированных точек будут высоко оценены.
справа 

Подробнее здесь: [url]https://stackoverflow.com/questions/79825133/how-can-i-make-my-solution-to-converge-for-determining-fixed-point-for-nonlinear[/url]