Мне нужно смоделировать движение стального стержня через магнит

# Сохраните как simulate_magnetization.py и запустите: python simulate_magnetization.py

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

mu0 = 4*np.pi*1e-7

params = {'Ms':1.6e6, 'a':1.0e5, 'k':5.0e4, 'c':0.2, 'alpha':1e-3}

sigma = 1.0e6

rod_thickness = 0.01

magnet_peak_B = 0.5

magnet_center = 0.0

magnet_sigma = 0.02

def B_external(x):

return magnet_peak_B np.exp(-0.5 ((x - magnet_center)/magnet_sigma)**2)

def M_anh(H_eff, Ms, a):

return Ms*np.tanh(H_eff/a)

tau = mu0 sigma (rod_thickness/2)**2

print("tau (s) =", tau)

def integrate_JA_time(H_ts, dt, params, tau):

Ms = params['Ms']; a = params['a']; k = params['k']; c = params['c']; alpha = params['alpha']

M = 0.0

M_hist = np.zeros_like(H_ts)

for i, H_ext in enumerate(H_ts):

dH_dt = 0.0 if i==0 else (H_ts[i]-H_ts[i-1])/dt

H_eddy = -tau * dH_dt

H_eff = H_ext + H_eddy

Man = M_anh(H_eff + alpha*M, Ms, a)

delta = np.sign(dH_dt) if dH_dt != 0 else 1.0

denom = (k delta - alpha (Man - M))

dM_dH = 0.0 if abs(denom) < 1e-12 else (Man - M) / denom

dM_dt = dM_dH * dH_dt

M = M + dM_dt * dt

M = c*Man + (1-c)*M

M_hist[i] = M

return M_hist

velocities = np.array([0.01,0.02,0.05,0.1,0.2,0.5])

results = []

for v in velocities:

x_start, x_end = -0.08, 0.08

total_time = abs(x_end-x_start)/v

nsteps = max(2000, int(total_time/1e-5))

t = np.linspace(0, total_time, nsteps)

dt = t[1]-t[0]

x_t = x_start + v*t

B_t = B_external(x_t)

H_t = B_t / mu0

M_t = integrate_JA_time(H_t, dt, params, tau)

results.append({'velocity_m_s':v, 'mean_abs_M_A_m':np.mean(np.abs(M_t)), 'peak_abs_M_A_m':np.max(np.abs(M_t))})

plt.figure()

plt.plot(t, M_t/params['Ms'])

plt.xlabel("time (s)")

plt.ylabel("M(t)/Ms")

plt.title(f"v={v} m/s")

plt.grid(True)

plt.show()

df = pd.DataFrame(results)

print(df)

df.to_csv("magnetization_vs_velocity_improved.csv", index=False)

print("Saved: magnetization_vs_velocity_improved.csv")