ARTICLE AD BOX
I'm trying to plot nullclines for a nonlinear system and determine fixed points from the same. However till now my code is not good enough to generate the plots. Presently I'm only able to get a blank plot only with various guess values tried for finding the root.
My objective is to generate the plots for g1 and I such that the intersecting points will be the fixed points whose coordinates are to be studied. Some coordinates as needed for the same are imported from a matlab file path included in the code.
Please advice me the corrections/fixes to be made. My code is as below
import sympy as sp from sympy import nsolve, Abs, symbols import numpy as np from scipy.optimize import brentq from pymatreader import read_mat import matplotlib.pyplot as plt import pandas as pd import traceback import math import warnings warnings.filterwarnings("ignore", category=RuntimeWarning) from matplotlib.lines import Line2D from sympy import Heaviside # --- Load MAT file (adjust path) --- file_path = r'E:\data&plot\data 3\set5\data1_x0_03.mat' data = read_mat(file_path) data1 = data['data1'] vars_ = data1['variable'] consts = data1['constant'] t = 9999 j = 99 i = 9999 # --- Pull constants ---- H0_val = consts.get('H0') 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') # note: your const key '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 = consts.get('critgradpressure', consts.get('critgradpressure', None)) Q_val = float(vars_['Q'][i, j]) I_val = float(vars_['turbulence'][i, j]) g1_val = float(vars_['gradpressure'][i, j]) #print(mu_c_val) # check for missing values missing = [k for k,v in [ ('chi0', chi0_val), ('chi1', chi1_val), ('D0', D0_val), ('D1', D1_val), ('alphachi', alpha_chi_val), ('alphaD', alpha_D_val), ('chi_growth/chi_', chi_val), ('beta_suppression/beta_', beta_val), ('lambda_suppress/lambda', lam_val), ('critgradpressure', mu_c_val) ] if v is None] if missing: raise RuntimeError("Missing constants in MAT file: " + ", ".join(missing)) # Convert to floats (some matlab values may already be floats) H0_val = float(H0_val) chi0_val = float(chi0_val) chi1_val = float(chi1_val) D0_val = float(D0_val) D1_val = float(D1_val) alpha_chi_val = float(alpha_chi_val) alpha_D_val = float(alpha_D_val) chi_val = float(chi_val) beta_val = float(beta_val) lam_val = float(lam_val) mu_c_val = float(mu_c_val) #print(mu_c_val) # Using v_eQ = -g1**2 => v_eQ**2 = g1**4 tot_vals = [] tot = [] I_null=[] g1_safe=[] def g1_dot_nullcline (g1_val, Q_val, chi0_val, chi1_val, alpha_chi_val, I_expr, chi_val, beta_val, lam_val, mu_c_val, H0_val, tot, tot_vals): g1_val = np.asarray(g1_val) print(g1_val) g1_safe = np.where(np.abs(g1_val) < 1e-12, 1e-12, g1_val) v_eQ_sq = ( -g1_safe**2 )**2 # (-g1^2)^2 = g1^4 A = Q_val*(chi0_val + chi1_val*I_val / (1 + alpha_chi_val*(-g1_safe**2)))**(-1) term1 = np.abs(-g1_safe - mu_c_val) term = term1 * (term1 > 0) # Heaviside Θ term_t = term - lam_val * v_eQ_sq B = chi_val*term_t-beta_val*I_val**2 tot_val=A*B-H0_val tot_vals.append(tot_val) return (tot_vals) for g1 in g1_safe: # Try solving in a reasonable I-range try: root = brentq(nullcline_eq, -30, 30, args=(g1,)) I_null.append(root) except: I_null.append(np.nan) I_null = np.array(I_null) # --- Plot --- plt.figure(figsize=(8,6)) plt.plot(g1_safe, I_null, lw=3) plt.xlabel("g₁") plt.ylabel("I (nullcline)") plt.title("g₁-nullcline: I(g₁) from g1_dot = 0") plt.grid(True) plt.show() print("total=", g1_dot_nullcline(g1_val, Q_val, chi0_val, chi1_val, alpha_chi_val, I_val, chi_val, beta_val, lam_val, mu_c_val, H0_val, tot, tot_vals)) g1_val = np.asarray(g1_val) print(tot_vals)