Хабр Курсы для всех
РЕКЛАМА
Большая витрина: от крупнейших школ до частных авторов. Сравнивайте по цене, длительности, формату и выбирайте самый подходящий курс!
Поток Научпоп доступен 24/7 благодаря поддержке друзей Хабра
Комментарии 65
Закреплённые комментарии
Общий подход в такого рода вещах примерно такой. Запускаете моделирование и долго на него смотрите. Спрашиваете себя, а понимаю ли я, что здесь происходит. Например, на картинке, подписанной "Волна на плоскости", никакой волны не видно. Область, в которой задано начальное отклонение, как-то себе крутится, а все остальное остается в покое. На первый взгляд, на границе, отделяющей область начального отклонения, лапласиан магнитного момента должен быть весьма значительным, что создает эффективное поле с компонентой нормальной исходному магнитному моменту. Значит, около области исходного отклонения магнитные моменты тоже должны начать крутиться и т.д. Т.е. в самом деле должно получиться распространяющееся возмущение, а его нет. В чем дело?
По картинкам, показанным в тексте, можно много вопросов позадавать, но начать нужно с самых базовых, самых простых. Скорость распространения возмущения, время жизни начального возмущения (да и можно ли поставить вопрос о времени жизни, если нет затухания?), насколько физично задавать начальное возмущение в том виде, как это сделано в модели, а как проверить, что модель записана правильно и нигде в коде никаких минусов не потерялось, и т.д., и т.п.
Ромбическая форма распространения волны как бы намекает что во всём этом скрыты артефакты вычислений и результатам не стоит верить на 100%.
Автор прямо так гиперактивно на меня ссылается, аж жуть.
Поэтому дисклеймер: его работы к моим никакого содержательного отношения не имеют!
А насколько мои анимации соответствуют действительности?
У вас ни масштабы времени не описаны, ни шаг решётки, ни остальные параметры — как они соотносятся с реальностью? Разберитесь, подставьте, сравните с реальностью, которую можно найти в статьях.
А то вы с темы на тему прыгаете, моделируете всё что под руку попадается, без глубокого погружения и анализа, но вместе с тем мечтаете в комментариях о публикациях. Любое исследование — это прежде всего систематическая долбёжка одной проблемы с тщательным и обоснованным выбором методов, а не накидывание красивых картинок спорной реалистичности, а не прыжки от распределения температур к ОТО, потом к полям магнитов, затем к гидродинамике, затем к уравнениям ЛЛГ. Как некоторые радиолюбители‑паяльщики, которым прежде всего важен процесс пайки — а сборка и настройка готовой конструкции откладывается в долгий ящик.
Определитесь сперва, что вам интересно, и потом уже вгрызайтесь в эту тему.
В коде чётко прописаны все ключевые параметры среды и сетки:
# ================= ПАРАМЕТРЫ (Хабр‑режим) =================
Nx, Ny = 64, 64 # сетка 64x64
dx = 10e-9 # шаг сетки 10 нм (крупнее → устойчивее и быстрее)
dt = 5e-13 # шаг по времени 0.5 пс
steps = 1200 # сколько шагов считать
save_every = 15 # сохранять каждый 15‑й кадр
gamma = 2.211e5 # гиромагнитное отношение
alpha = 0.005 # небольшое затухание (чтобы волна не «звенела» вечно)
A_ex = 1.3e-11 # обменная константа
Ms = 8e5
mu0 = 4 np.pi 1e-7
H_ext = 4.0e4 # УВЕЛИЧЕННОЕ поле (4 Т) — чтобы динамика была видна
Ну а хоть на какой-нибудь реальный магнетик эти числа похожи?
Эти параметры очень близко описывают пермаллой (Permalloy) — сплав на основе железа и никеля (обычно ~80 % Ni, 20 % Fe).
Ms=8⋅105 А/м — попадает в диапазон 7.5–8.5⋅105 А/мЭто классика для Ni80Fe20.
Aex=1.3⋅10−11 Дж/мA— тоже типично для пермаллоя (часто берут 1.3 или 1.4⋅10−11.
γ=2.211⋅105 рад/(с ⋅А/м)— соответствует гиромагнитному отношению для намагниченности насыщения пермаллоя.
α=0.005 — реалистичное значение для качественного пермаллоя в тонких плёнках. В образцах с хорошей кристаллической структурой и малой дефектностью затухание бывает порядка 0.003–0.008
Эту статью я специально написал под тему ваших научных интересов
Что касается разнообразия тем для моих статей, то в одном из постов вы писали следующее:
Поскольку физик‑теоретик — понятие очень абстрактное, размытое и растяжимое, заниматься по научной тематике доводилось различными задачами. Гидродинамика, нелинейные процессы, спектральный и вейвлет‑анализ, немного биомедицинских работ. Это всё немного пылится в прошлом, а сейчас основной темой стали спиновая динамика и всяческие волны намагниченности в твёрдых телах.
Источник:https://habr.com/ru/articles/737624/
Эту статью я специально написал под тему ваших научных интересов: вы же недавно защитили докторскую диссертацию по этой тематике (спиновая динамика и волны намагниченности). Материал интересный и уникальный: до меня такого никто на Хабре не делал. Параметры среды соответствуют реальному веществу- пермаллою, мне было очень интересно их двигать и смотреть на полученные анимации (как меняется динамика). Думаю, вам тоже интересно их посмотреть...
Сейчас я стараюсь делать качественный научно-популярный контент с полноценными анимациями, кодом, формулами и списком литературы. Кстати, нейросеть (Алиса) мне в этом отлично помогает: она мне пишет код по теме, которую я попросил, остаётся только накидать текста и формул с гиперссылками- вуаля, и статья готова за всего 1-3 дня.
Современные технологии шагнули вперёд. То, чему люди раньше учились несколько лет в ВУЗах, сейчас с нейросетями спокойно может делать толковый человек даже без специального образования...
Нейросети не заменят понимания
Ну я же примерно понимаю моделируемые физические явления (сдал ЕГЭ по физике на 97 баллов).
Я же написал :толковый человек.
Сданный егэ показывает лишь, что вы имеете достаточную базу для продолжения профильного обучения.
Вы же просто хватаете очередные попавшиеся вам под руку уравнения, и пишете для них код, который еле работает, не зная даже элементарных особенностей схем и оценки устойчивости, и прикрываясь числом Куранта, которое, конечно, критерий устойчивости, но оно ничего не говорит о монотонности и консервативности схем. Вы же даже не в состоянии самостоятельно оценить корректность результатов, раз постоянно спрашиваете об этом. Результат нужно защитить, обосновать, а не жалостливо просить посмотреть на него.
Возьмите десятитомник Ландау и посмотрите там задачи после параграфов, задачи к курсу лекций Фейнмана, задачник по теорфизике Гречко или Белоусова, сборник качественных задач «Понимаете ли вы физику» Капицы. Если сможете хотя бы в 10–20 процентах задач по какому‑нибудь разделу физики сделать первый шаг без генератора случайных слов нейросети — да, вы что‑то понимаете в этом разделе.
Преодолейте этот уровень, поймите, насколько иронично выглядят ваши расчёты во всех статьях — если получится, можно продолжать. И найти какого‑нибудь научрука, который разделит ваш энтузиазм.
А как можно обосновать, что при постепенно увеличении внешнего поля колебания сначала ускорились потом значительно замедлились и стали хаотичными, а потом опять ускорились и стали синхронными? Напишите пожалуйста, я вставлю в статью.
Вы автор результата. Вам обосновывать происходящее, а также доказывать, что этого никто раньше не получал
При постепенном увеличении амплитуды внешнего поля в системе, описываемой уравнением ЛЛГ, наблюдается последовательность качественно различных режимов колебаний. На начальном этапе слабое поле действует как малая модификация эффективного поля, задаваемого материалом (включая анизотропию и размагничивающие поля). В этих условиях система демонстрирует ускорение колебаний: внешнее поле начинает резонансно «подхватывать» собственные моды прецессии намагниченности, что приводит к росту амплитуды и частоты отклика.
По мере дальнейшего роста поля нелинейные эффекты выходят на первый план. Эффективный потенциал, определяющий динамику вектора намагниченности, усложняется: возникают дополнительные минимумы, барьеры, области с сильной зависимостью от угла между намагниченностью и полем. В такой ситуации даже небольшое возмущение может привести к тому, что траектории в фазовом пространстве начинают пересекать границы различных резонансных областей. Это порождает сложную интерференцию мод, бифуркации и, как следствие, переход к хаотическому режиму: колебания становятся нерегулярными, их частота и амплитуда перестают быть предсказуемыми. academia.edu +1
При ещё более сильных полях система может выйти на новый упорядоченный режим. В определённых областях параметра (например, когда внешнее поле доминирует и задаёт чёткий «ведущий» ритм) снова возникает условие для устойчивой синхронизации. Система «захватывается» в новый резонансный режим: колебания стабилизируются, их частота и фаза выстраиваются в регулярный, синхронный ритм. Важно отметить, что эта синхронная стадия не обязательно возвращает систему к исходной простой картине — скорее, возникает новый, возможно, более сложный, но упорядоченный режим.» pps.kaznu.kz +1
А если без ИИ, а своей головой?
Вначале внешнее поле невелико, доминирующим вкладом в H_eff становится само внешнее поле. В этом случае можно приближённо считать, что H_eff ≈ μ₀H_ext (с учётом магнитной постоянной). Тогда прецессия вектора M происходит с частотой, близкой к ларморовской, и частота колебаний линейно от него зависит.
При увеличении поля вклад других составляющих H_eff (анизотропия, обмен) становится сопоставимым с вкладом внешнего поля. Система становится чувствительной к начальным условиям и малым возмущениям. Решение уравнения ЛЛГ в этой области перестаёт быть простым гармоническим колебанием — появляются сложные, нелинейные траектории.
При достаточно большом внешнем поле система может выйти на новый устойчивый режим.
Я добавил ещё оценку устойчивости при помощи показателя Ляпунова, чтобы обосновать это всё:
Устойчивость динамики намагниченности оценивалась по локальному показателю Ляпунова, рассчитанному численно путём сравнения двух близких траекторий в фазовом пространстве. При малых значениях внешнего поля показатель λ<0, что соответствует устойчивым регулярным колебаниям (режим ускорения). В промежуточной области наблюдается переход к λ>0, указывающий на возникновение хаотической динамики и чувствительность к начальным условиям. При высоких полях показатель снова становится отрицательным (λ<0), что свидетельствует о восстановлении устойчивости и переходе к режиму вынужденной синхронизации.
При этом показатель Ляпунова на всём диапазоне колебался около нуля, но на подавляющем большинстве интервалов оставался отрицательным.
Теперь эту работу можно защитить перед вами, известным специалистом по спиновой динамике и доктором физико-математических наук. Задавайте вопросы.
А сейчас можно ли считать полученный результат полностью обоснованным и защитить его перед вами??
Скажите пожалуйста, в ходе научных исследований по спиновой динамике делали ли вы что-то подобное (писали код, делали такую же симуляцию)?
Я в этом отношении классический теоретик. Первопринципные модели и аналитические результаты в первую очередь. Численный счёт сильно потом.
Ну и наконец, могут ли уравнения ЛЛГ, Навье-Стокса, теории относительности, температурного и магнитного поля и т. д. вдруг выпасть на школьном ЕГЭ по физике?
Вот я сделал процедуру обезразмеривания и результат теперь другой:
Теперь всё реалистично
Общий подход в такого рода вещах примерно такой. Запускаете моделирование и долго на него смотрите. Спрашиваете себя, а понимаю ли я, что здесь происходит. Например, на картинке, подписанной "Волна на плоскости", никакой волны не видно. Область, в которой задано начальное отклонение, как-то себе крутится, а все остальное остается в покое. На первый взгляд, на границе, отделяющей область начального отклонения, лапласиан магнитного момента должен быть весьма значительным, что создает эффективное поле с компонентой нормальной исходному магнитному моменту. Значит, около области исходного отклонения магнитные моменты тоже должны начать крутиться и т.д. Т.е. в самом деле должно получиться распространяющееся возмущение, а его нет. В чем дело?
По картинкам, показанным в тексте, можно много вопросов позадавать, но начать нужно с самых базовых, самых простых. Скорость распространения возмущения, время жизни начального возмущения (да и можно ли поставить вопрос о времени жизни, если нет затухания?), насколько физично задавать начальное возмущение в том виде, как это сделано в модели, а как проверить, что модель записана правильно и нигде в коде никаких минусов не потерялось, и т.д., и т.п.
Обосновать полученный результат можно так:
Вначале внешнее поле невелико, доминирующим вкладом в H_eff становится само внешнее поле. В этом случае можно приближённо считать, что H_eff ≈ μ₀H_ext (с учётом магнитной постоянной). Тогда прецессия вектора M происходит с частотой, близкой к ларморовской, и частота колебаний линейно от него зависит.
При увеличении поля вклад других составляющих H_eff (анизотропия, обмен) становится сопоставимым с вкладом внешнего поля. Система становится чувствительной к начальным условиям и малым возмущениям. Решение уравнения ЛЛГ в этой области перестаёт быть простым гармоническим колебанием — появляются сложные, нелинейные траектории.
При достаточно большом внешнем поле система может выйти на новый устойчивый режим.
Вначале внешнее поле невелико, доминирующим вкладом в H_eff становится само внешнее поле.
А это действительно так? Все же, на первый взгляд, зависимость (скажем, вертикальной, по картинке, компоненты) намагниченности от координаты выглядит разрывной функцией: вот здесь, за пределами области начального отклонения, она ноль, а вот здесь - совсем не ноль. Уже первая производная такой функции равна бесконечности, а про вторую и говорить не надо. Т.е. в пространственной области перехода от возмущенной намагниченности к исходной эффективное поле может быть как угодно большим. Более точно, оно будет обратно пропорционально квадрату ширины этой области. В используемой модели минимальный масштаб - шаг решетки, т.е. 10 нм. Т.к. моделирование не позволяет увидеть эффект разрыва намагниченности, можно сделать вывод, что шаг слишком большой для того, чтобы вообще увидеть эффект обменного взаимодействия: взаимная ориентация намагниченности в соседних узлах расчетной решетки может быть какой угодно. Или же какая-то другая неполадка происходит в коде.
Как, вообще, отвечать на вопрос о выборе шага решетки? Есть правило "большого пальца" - уменьшили шаг решетки вдвое, посмотрели на результат. Если серьезно не поменялся, то все в порядке. Проблема с этим подходом в том, что он работает только когда шаг с самого начала был примерно таким, каким нужно. А когда он совсем не в тему, то уменьшение в два раза приведет только лишь к тому, что он станет чуть меньше не в тему. Есть у меня связанная с этим история, которая продолжалась несколько месяцев, но это для другого раза.
Мы пришли к одному очень-очень важному выводу. Моделирование должно характеризоваться не материальными, а физическими параметрами. Получается некоторый парадокс: мы хотим что-то промоделировать, но для того, чтобы это сделать, нам уже нужно знать результаты моделирования. И это действительно так. За всем этим стоят свои большие теории, но, к счастью, туда пока можно не соваться.
Первый шаг, который делают в направлении определения физических параметров, максимально обезразмеривают уравнение. Например, намагниченность представляется как единичный вектор, помноженный на некоторую величину, время - безразмерный параметр, умноженный на характерный масштаб, то же с координатой. Записывают уравнение в новых координатах и выбирают характерные масштабы так, чтобы итоговое уравнение имело как можно меньше свободных параметров. Процедура несложная, но проводить ее лучше в явном виде, особенно если уравнение нелинейное, как уравнене Ландау-Лифшица. В модуле, который занимается собственно вычислением динамики, вот этого блока
gamma = 2.211e5 # гиромагнитное отношение
alpha = 0.005 # небольшое затухание (чтобы волна не «звенела» вечно)
A_ex = 1.3e-11 # обменная константа
Ms = 8e5
mu0 = 4 * np.pi * 1e-7быть не должно. Все эти константы в уравнении, представленном для вычислений, должны быть единицы, а нанометры с наносекундами должны получаться только на стадии вывода результатов для человека.
Тем не менее в коде, представленном в заметке, этот блок присутствует (это плохой знак). Проверяем насколько последовательно материальные параметры представлены в вычислительной части. У меня глаз сразу зацепился за
# Нормировка
norm = np.sqrt(np.sum(M**2, axis=2, keepdims=True))
norm[norm == 0] = 1.0
M = M / normТ.е. намагниченность фактически считается единичным вектором, тогда как исходное уравнение записано для физической намагниченности. Поскольку уравнение Ландау-Лифшица нелинейно, на такого рода расхождения нужно смотреть особенно пристально. В данном случае создается впечатление, что эффективно обменное взаимодействие почти на шесть порядков понизилось. Поэтому вполне можно ожидать, что для такого эффективного взаимодействия (т.е. для гипотетического материала) шаг решетки оказался гигантским, и спиновые волны он может не разглядеть.
А что вы можете сказать про оценку устойчивости по Ляпунову? Всё верно?
При чём тут устойчивость по Ляпунову, если вы с устойчивостью счёта разобраться не можете?
Толковые люди идут по толковому пути. Вы утверждаете, что это студенческая работа, в профиле пишете, что стажёр‑энергетик. Если вы правда вхожи в какой‑нибудь университет, то наверняка у вас есть доступ к какой‑нибудь кафедре математики и кафедре физики.
Вы всё время задаёте вопросы, которые сперва должен оценивать научный руководитель, и попутно должен учить вас отвечать на эти вопросы самостоятельно. В этом и есть одна из ключевых сторон квалификации исследователя. Опять же — понимание происходящего в изучаемом явлении. Мало назваться «толковым человеком».
Ничего не могу сказать. Вероятно, формулы посчитаны верно, но на какой вопрос они дают ответ, сказать непросто. Конечно, слова какие-то наговорить можно, но их ценность невелика.
У меня к такого рода ситуациям подход физика. Если я не вижу какого-то элементарного ожидаемого эффекта (в данном случае, спиновых волн), то любые методики, которые говорят, что все в порядке, получают характеристику "не помогают разобраться в происходящем".
Ваши численные результаты говорят, что (в относительно низких внешних полях) в пермаллое нет спиновых волн, у уравнения Ландау-Лифшица-Гильберта нет волновых решений. Это один из тех экстраординарных результатов, что требуют экстраординарных доказательств.
Приведу аналогию. Вам показывают моделирование двух блоков, связанных пружиной, без всяких внешних сил. Один блок колеблется, другой стоит как влитой. Ваша реакция? А потом вам говорят, что вот, пожалуйста, устойчивость по Ляпунову. Ваша реакция в этом случае?
Что экстраординарного в этом результате?
В относительно низких полях колебания все равно есть, а значит и волны тоже. Возможно, они распространяются очень медленно и мы их на гифке не успеваем увидеть.
А как можно отличить ситуации "волны есть, но мы их не видим" и "волн нет, и мы их не видим"?
Переформулирую экстраординарный результат "уравнение Ландау-Лифшица с обменным взаимодействием допускает решения с пространственно разрывной намагниченностью". Такое решение явно показано на картинке. Экстраординарность в том, что само представление об обменном взаимодействии в уравнении Ландау-Лифшица получается в предположении о непрерывности.
Другой экстраординарный результат, связанный с разрывностью, более тонкий. Уравнение Ландау-Лифшица допускает топологические солитоны - доменные стенки. В строгом отсутствии анизотропии, но при наличии внешнего поля - это 360-градусные доменные стенки. Используя ваш код, предположительно можно "продемонстрировать", что начальное состояние с такой стенкой эволюционирует в состояние без нее. Достаточно задать такой разворот, включить затухание побольше (наличие затухание не влияет на существование топологических солитонов) и все магнитные моменты по скучным траекториям свалятся к внешнему полю. Правда, возможен артефакт соизмеримости начального условия и решетки, где выживет слой с намагниченностью, направленной против поля.
Решите проблему с разрывностью намагниченности.
На какой конкретно картинке видно пространственно разрывную намагниченность? Я не вижу.
На картинке "Волна на плоскости", речь все еще про нее. Нарисуйте, скажем, зависимость у-компоненты намагниченности от у при х = 0. Если не убеждает, нарисуйте ее вторую производную по y (это дает вклад в эффективное поле).
Если можно взять вторую производную(как вы утверждаете), то намагниченность непрерывна по y.
Речь идет о работе с массивом данных, который получился у вас в программе. Т.е. с массивом, который изображен на картинке "Волна на плоскости". У вас в программе от него Лапласиан вычисляется. Оставить от него только вторую производную по у не должно составить труда.
Всё непрерывно, никакой пространственно разрывной намагниченности нет.
Это, очевидно, никак не соответствует ни картинке "Волна на плоскости", ни коду, который приведен в заметке.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# --- Параметры пермаллоя (Ni80Fe20) ---
Ms = 800e3 # А/м, намагниченность насыщения
A_ex = 13e-12 # Дж/м, обменная константа
alpha = 0.005 # безразмерный параметр затухания
gamma = 2.211e5 # рад/(с·А/м), гиромагнитное отношение
# --- Сетка (одномерная вдоль y) ---
L = 200e-9 # длина образца (200 нм)
N = 401 # нечётное число точек (центр ровно посередине)
y = np.linspace(-L/2, L/2, N)
dy = y[1] - y[0]
# --- Внешнее поле (вдоль x) ---
H_ext_x = 4e4 # А/м (поле, которое сдвигает доменную стенку)
# --- Начальное условие: доменная стенка (Блоховская) ---
# Толщина стенки без анизотропии: delta ~ sqrt(A/K), K=0 → берём условную толщину
delta = 20e-9 # 20 нм — реалистичная толщина стенки для пермаллоя
theta_init = np.pi/2 * (1.0 + np.tanh(y / delta))
mx = np.sin(theta_init)
my = np.cos(theta_init)
mz = np.zeros_like(y)
M = np.stack([mx, my, mz], axis=0) # M[0]=Mx, M[1]=My, M[2]=Mz
# --- Эффективное поле: обмен + внешнее ---
def compute_Heff(M_norm):
mx, my, mz = M_norm
# Вторая производная через np.gradient (дважды)
d2mx_dy2 = np.gradient(np.gradient(mx, dy), dy)
d2my_dy2 = np.gradient(np.gradient(my, dy), dy)
d2mz_dy2 = np.gradient(np.gradient(mz, dy), dy)
pref = (2.0 * A_ex) / (Ms**2)
H_exchange_x = pref * d2mx_dy2
H_exchange_y = pref * d2my_dy2
H_exchange_z = pref * d2mz_dy2
H_ext = np.array([H_ext_x, 0.0, 0.0])
Heff = np.stack([H_exchange_x, H_exchange_y, H_exchange_z], axis=0) + H_ext[:, None]
return Heff
# --- Правая часть уравнения ЛЛГ (Gilbert form) ---
def llg_rhs(M_norm, t):
mx, my, mz = M_norm
Heff = compute_Heff(M_norm)
Hx, Hy, Hz = Heff
# M × H_eff
cx = my * Hz - mz * Hy
cy = mz * Hx - mx * Hz
cz = mx * Hy - my * Hx
# M × (M × H_eff)
c2x = my * cz - mz * cy
c2y = mz * cx - mx * cz
c2z = mx * cy - my * cx
dMdt_x = -gamma * cx + alpha * gamma * c2x
dMdt_y = -gamma * cy + alpha * gamma * c2y
dMdt_z = -gamma * cz + alpha * gamma * c2z
return np.stack([dMdt_x, dMdt_y, dMdt_z], axis=0)
# --- Интегрирование (явный Эйлер с нормировкой) ---
dt = 5e-14 # 50 фс — маленький шаг для устойчивости
t_max = 2e-9 # 2 нс
save_every = 4 # сохранять каждый 4-й кадр
n_steps = int(t_max / dt)
times = []
My_profiles = []
M_curr = M.copy()
t = 0.0
for step in range(n_steps):
dM = llg_rhs(M_curr, t)
M_curr += dM * dt
# Нормировка |M| = 1 (строго, без NaN)
norm = np.sqrt(M_curr[0]**2 + M_curr[1]**2 + M_curr[2]**2)
# Защита от нулевого вектора (на всякий случай)
norm[norm == 0] = 1.0
M_curr /= norm
t += dt
if step % save_every == 0:
times.append(t)
My_profiles.append(M_curr[1].copy())
times = np.array(times)
My_profiles = np.array(My_profiles)
# --- Анимация ---
fig, ax = plt.subplots(figsize=(8, 5))
ax.set_xlim(-L/2 * 1e9, L/2 * 1e9)
ax.set_ylim(-1.1, 1.1)
ax.set_xlabel('y (нм)')
ax.set_ylabel(r'$M_y / M_s$')
ax.set_title('Динамика доменной стенки в пермаллое (x=0)')
ax.grid(True, alpha=0.3)
line, = ax.plot([], [], 'b-', linewidth=2)
time_text = ax.text(0.05, 0.95, '', transform=ax.transAxes, fontsize=12, verticalalignment='top')
def init():
line.set_data([], [])
time_text.set_text('')
return line, time_text
def animate(i):
y_nm = y * 1e9
line.set_data(y_nm, My_profiles[i])
time_text.set_text(f't = {times[i]*1e9:.2f} нс')
return line, time_text
ani = FuncAnimation(fig, animate, frames=len(times), init_func=init, blit=True, interval=30)
plt.show()
# Чтобы сохранить в GIF (если установлен pillow):
ani.save('permalloy_domain_wall.gif', writer='pillow', fps=15)
Вот код как был получен этот график.
Я предполагал что-то такое. Главная ошибка исходного кода перенесена в этот. Началась она с методической: не обезразмерено решаемое уравнение. Перешела в вычислительную: неправильно вычисляется обменное взаимодействие.
Повторюсь. В уравнении Ландау-Лифшица присутствует физическая намагниченность, т.е. вектор с длиной 
. В вашем коде намагниченность представлена единичным вектором 
. По меньшей мере, в части, которая вычисляет эффективное поле, должно быть не
pref = (2.0 * A_ex) / (Ms**2)а
pref = (2.0 * A_ex) / MsВ правильном же коде должно быть вообще что-то вроде
pref = 1В хорошем коде можно ожидать свои хитрости, но о них вести речь сейчас бессмысленно.
Для новой постановки задачи физические несоответствия увидеть сложнее. Если убрать затухание, то будет проще увидеть, что динамика сохраняет структуру исходной доменной стенки, а она не должна этого делать. А если есть затухание, то нужно смотреть на 360-градусную доменную стенку. Тогда динамика должна произвольное такое начальное распределение приводить в "стандартную" форму. В вашем же коде, если подождать, возможно полезут разрывы. Хотя тут возможны варианты, т.к. в новом коде шаг решетки в 20 раз более мелкий. Навскидку сказать достаточно этого или нет моего внутреннего калькулятора не хватает. Вроде бы все еще шаг слишком здоровый. В правильно обезразмеренном уравнении все это было бы очевидным.
Вам же все равно LLM код пишет? Попросите ее обезразмерить (в итоговом коде внутри не должно быть материальных констант вообще) и написать соответствующие преобразующие функции.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation
# =========================================================
# ПАРАМЕТРЫ ПЕРМАЛЛОЯ (размерные)
# =========================================================
Ms = 800e3
A_ex = 13e-12
alpha = 0.005
gamma = 2.211e5
# =========================================================
# ГЕОМЕТРИЯ И СЕТКА
# =========================================================
L = 200e-9
N = 401
y = np.linspace(-L/2, L/2, N)
dy = y[1] - y[0]
# Безразмерные единицы
l0 = 10e-9 # единица длины (10 нм)
t0 = 1.0 / (gamma * Ms) # единица времени (~0.57 пс)
y_dimless = y / l0
dy_dimless = dy / l0
# Безразмерный обменный коэффициент
J_dimless = (2.0 * A_ex) / (Ms**2 * l0**2) * t0 * gamma
# Поле размагничивания и внешнее поле (безразмерные)
N_perp = 0.25
H_ext_x_dimless = 4e4 / Ms
# =========================================================
# НАЧАЛЬНОЕ УСЛОВИЕ
# =========================================================
center_idx = N // 2
width_dimless = 4.0
mx_init = np.tanh((y_dimless[center_idx] - y_dimless) / (width_dimless / 2.0))
my_init = np.zeros_like(y_dimless)
mz_init = np.zeros_like(y_dimless)
norm = np.sqrt(mx_init**2 + my_init**2 + mz_init**2)
norm[norm == 0] = 1.0
mx_init /= norm
my_init /= norm
mz_init /= norm
m = np.stack([mx_init, my_init, mz_init], axis=0) # [3, N], |m|=1
# =========================================================
# ЭФФЕКТИВНОЕ ПОЛЕ (БЕЗРАЗМЕРНОЕ)
# =========================================================
def compute_h_eff(m_vec):
mx, my, mz = m_vec
d2mx = np.gradient(np.gradient(mx, dy_dimless), dy_dimless)
d2my = np.gradient(np.gradient(my, dy_dimless), dy_dimless)
d2mz = np.gradient(np.gradient(mz, dy_dimless), dy_dimless)
h_exchange_x = J_dimless * d2mx
h_exchange_y = J_dimless * d2my
h_exchange_z = J_dimless * d2mz
h_demag_y = -N_perp * my
h_ext = np.array([H_ext_x_dimless, 0.0, 0.0])
h_eff = np.stack(
[h_exchange_x + h_ext[0], h_exchange_y + h_demag_y, h_exchange_z + h_ext[2]],
axis=0
)
return h_eff
# =========================================================
# ПРАВАЯ ЧАСТЬ УРАВНЕНИЯ ЛЛГ (БЕЗРАЗМЕРНАЯ)
# =========================================================
def llg_rhs(m_vec, tau):
mx, my, mz = m_vec
h = compute_h_eff(m_vec)
hx, hy, hz = h
cx = my * hz - mz * hy
cy = mz * hx - mx * hz
cz = mx * hy - my * hx
c2x = my * cz - mz * cy
c2y = mz * cx - mx * cz
c2z = mx * cy - my * cx
dm_dt_x = -cx - alpha * c2x
dm_dt_y = -cy - alpha * c2y
dm_dt_z = -cz - alpha * c2z
return np.stack([dm_dt_x, dm_dt_y, dm_dt_z], axis=0)
# =========================================================
# ИНТЕГРИРОВАНИЕ (РЕЛАКСАЦИЯ + ДВИЖЕНИЕ)
# =========================================================
dt_tau = 0.05
# Релаксация
tau_relax = 50.0
n_steps_relax = int(tau_relax / dt_tau)
m_curr = m.copy()
for step in range(n_steps_relax):
dm = llg_rhs(m_curr, 0.0)
m_curr += dm * dt_tau
# Нормировки нет: структура уравнения сохраняет длину
# Движение
tau_move = 60.0
n_steps_move = int(tau_move / dt_tau)
save_every = 3
taus = []
my_profiles = []
mx_profiles = [] # для стрелки
for step in range(n_steps_move):
dm = llg_rhs(m_curr, step * dt_tau)
m_curr += dm * dt_tau
if step % save_every == 0:
taus.append(step * dt_tau)
my_profiles.append(m_curr[1].copy())
mx_profiles.append(m_curr[0].copy())
taus = np.array(taus)
my_profiles = np.array(my_profiles)
mx_profiles = np.array(mx_profiles)
times_real = taus * t0 # секунды
print(f"✅ Рассчитано кадров анимации: {len(taus)}")
print(f"[РЕЛАКСАЦИЯ] m_y в центре = {m_curr[1][N//2]:.3f} (не ноль — реалистично)")
# =========================================================
# ПОДГОТОВКА АНИМАЦИИ
# =========================================================
fig, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], linewidth=2, color='#1f77b4')
arrow = ax.arrow(0, 0, 0, 0, width=0.015, color='red', head_width=0.06, head_length=0.12)
ax.set_xlim(y[0]*1e9, y[-1]*1e9)
ax.set_ylim(-1.1, 1.1)
ax.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.4)
ax.set_xlabel('Координата y (нм)', fontsize=12)
ax.set_ylabel(r'$m_y$', fontsize=12)
title_text = ax.set_title('', fontsize=14)
ax.grid(True, linestyle=':', alpha=0.3)
def init():
line.set_data([], [])
arrow.remove()
global arrow
arrow = ax.arrow(0, 0, 0, 0, width=0.015, color='red', head_width=0.06, head_length=0.12)
title_text.set_text('')
return line, arrow, title_text
def update(frame):
y_nm = y * 1e9
my = my_profiles[frame]
mx = mx_profiles[frame]
line.set_data(y_nm, my)
# Удаляем старую стрелку и рисуем новую
arrow.remove()
center_y_nm = 0
center_my = my[N//2]
center_mx = mx[N//2]
# Стрелка: рисуем от центра немного вниз, чтобы не перекрывать график
arrow_y_start = -0.3
arrow_dx = center_mx * 20 # масштабируем стрелку
arrow_dy = center_my * 20
global arrow
arrow = ax.arrow(center_y_nm, arrow_y_start, arrow_dx, arrow_dy,
width=0.012, color='red', head_width=0.05, head_length=0.09)
t_ns = times_real[frame] * 1e9
title_text.set_text(f'Динамика доменной стенки: t = {t_ns:.2f} нс')
return line, arrow, title_text
ani = FuncAnimation(fig, update, frames=len(taus), init_func=init, blit=True, interval=40)
# Сохранение в GIF
gif_file = 'domain_wall_animation.gif'
ani.save(gif_file, writer='pillow', fps=20, dpi=200)
print(f"✅ Анимация GIF сохранена как {gif_file}")
# Сохранение в MP4 (нужен ffmpeg в PATH)
mp4_file = 'domain_wall_animation.mp4'
try:
ani.save(mp4_file, writer='ffmpeg', fps=20, dpi=200)
print(f"✅ Анимация MP4 сохранена как {mp4_file}")
except RuntimeError:
print("⚠️ Не удалось сохранить MP4: проверьте, установлен ли ffmpeg и добавлен ли в PATH.")
plt.close(fig)
Вот сделал обезразмеривание и убрал нормировку как вы просили. Результат такой же:
Все непрерывно.
Я выше описывал как обезразмеривание проводится. В новых обозначениях у вас должно быть
J_dimless = 1С временем процедура проведена верно, обратите внимание, что в llg_rhs уже нет gamma, такое же преобразование надо проделать и с пространственной координатой. В вычислении обменного взаимодействия должен быть чистый лапласиан. Т.е. должно быть не случайно выбранное l0 = 10e-9, а что-то осмысленное как в t0 = 1.0 / (gamma * Ms).
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io
# =========================================================
# 1. РАЗМЕРНЫЕ ПАРАМЕТРЫ (только для расчёта масштабов)
# =========================================================
Ms = 800e3
A_ex = 13e-12
alpha = 0.005
gamma = 2.211e5
# =========================================================
# 2. МАСШТАБЫ (без mu0 —consistent с исходным кодом)
# l0 = sqrt(2A / Ms^2) → J = 1
# t0 = 1 / (gamma * Ms)
# =========================================================
l0 = np.sqrt(2.0 * A_ex / Ms**2) # ≈ 6.374 нм
t0 = 1.0 / (gamma * Ms) # ≈ 5.654 пс
# =========================================================
# 3. БЕЗРАЗМЕРНЫЕ ВЕЛИЧИНЫ
# =========================================================
J = 1.0
h_ext_x = 4e4 / Ms # ≈ 0.05
L = 200e-9 / l0 # ≈ 31.39
N = 401
y = np.linspace(-L / 2, L / 2, N)
dy = y[1] - y[0]
# Начальное условие: доменная стенка
delta = 20e-9 / l0 # ≈ 3.14
theta = np.pi / 2 * (1.0 + np.tanh(y / delta))
mx = np.sin(theta)
my = np.cos(theta)
mz = np.zeros_like(y)
m = np.stack([mx, my, mz], axis=0)
# =========================================================
# 4. БЕЗРАЗМЕРНОЕ ЭФФЕКТИВНОЕ ПОЛЕ (J = 1)
# =========================================================
def compute_h(m_vec):
mx, my, mz = m_vec
d2mx = np.gradient(np.gradient(mx, dy), dy)
d2my = np.gradient(np.gradient(my, dy), dy)
d2mz = np.gradient(np.gradient(mz, dy), dy)
h = np.stack([d2mx, d2my, d2mz], axis=0)
h[0] += h_ext_x
return h
# =========================================================
# 5. БЕЗРАЗМЕРНОЕ УРАВНЕНИЕ ЛЛГ
# =========================================================
def llg_rhs(m_vec):
h = compute_h(m_vec)
mx, my, mz = m_vec
hx, hy, hz = h
cx = my * hz - mz * hy
cy = mz * hx - mx * hz
cz = mx * hy - my * hx
c2x = my * cz - mz * cy
c2y = mz * cx - mx * cz
c2z = mx * cy - my * cx
return np.stack([
-cx + alpha * c2x,
-cy + alpha * c2y,
-cz + alpha * c2z
], axis=0)
# =========================================================
# 6. ИНТЕГРИРОВАНИЕ: метод средней точки (RK2)
# k1 = f(m^n)
# k2 = f(m^n + dt/2 * k1)
# m^{n+1} = m^n + dt * k2
#
# Дрейф |m| за шаг — O(dt^4), за всё время — доли процента.
# Нормировка не нужна.
# =========================================================
dt = 0.01
tau_max = 2e-9 / t0 # ≈ 353.8
n_steps = int(tau_max / dt)
save_every = max(1, n_steps // 60)
taus, my_profiles, mx_profiles = [], [], []
m_curr = m.copy()
for step in range(n_steps):
k1 = llg_rhs(m_curr)
k2 = llg_rhs(m_curr + 0.5 * dt * k1)
m_curr += dt * k2
if np.any(~np.isfinite(m_curr)):
print(f"⚠️ NaN на шаге {step} — уменьшите dt")
break
if step % save_every == 0 or step == n_steps - 1:
taus.append(step * dt)
my_profiles.append(m_curr[1].copy())
mx_profiles.append(m_curr[0].copy())
taus = np.array(taus)
my_profiles = np.array(my_profiles)
mx_profiles = np.array(mx_profiles)
times_real = taus * t0
# Проверка нормы
norm_check = np.sqrt(m_curr[0]**2 + m_curr[1]**2 + m_curr[2]**2)
print(f"l0 = {l0*1e9:.3f} нм, t0 = {t0*1e12:.3f} пс, J = {J}")
print(f"Шагов: {n_steps}, кадров: {len(taus)}, время до {times_real[-1]*1e9:.2f} нс")
print(f"|m|: min={norm_check.min():.6f}, max={norm_check.max():.6f} (без нормировки)")
print(f"m_y в центре = {my_profiles[-1][N//2]:.3f}")
# =========================================================
# 7. АНИМАЦИЯ (GIF через Pillow)
# =========================================================
fig, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], linewidth=2, color='#1f77b4')
arrow_holder = [None]
y_nm = y * l0 * 1e9
ax.set_xlim(y_nm[0], y_nm[-1])
ax.set_ylim(-1.1, 1.1)
ax.axhline(0, color='k', linestyle='--', linewidth=1, alpha=0.4)
ax.set_xlabel('y (нм)', fontsize=12)
ax.set_ylabel(r'$m_y$', fontsize=12)
title = ax.set_title('', fontsize=14)
ax.grid(True, linestyle=':', alpha=0.3)
frames = []
for i in range(len(taus)):
line.set_data(y_nm, my_profiles[i])
if arrow_holder[0] is not None:
arrow_holder[0].remove()
cmx = mx_profiles[i][N // 2]
cmy = my_profiles[i][N // 2]
arrow_holder[0] = ax.arrow(
0, -0.3, cmx * 20, cmy * 20,
width=0.012, color='red',
head_width=0.05, head_length=0.09
)
title.set_text(f't = {times_real[i] * 1e9:.2f} нс')
buf = io.BytesIO()
fig.savefig(buf, format='png', dpi=100)
buf.seek(0)
frames.append(Image.open(buf).convert('RGB'))
plt.close(fig)
gif_file = 'domain_wall_nondim_RK2.gif'
if frames:
frames[0].save(
gif_file,
save_all=True,
append_images=frames[1:],
duration=50,
loop=0
)
print(f"✅ GIF сохранён: {gif_file}")Теперь всё полностью корректно:
результат-
стало быстрее, но всё непрерывно
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io
# =========================================================
# 1. РАЗМЕРНЫЕ ПАРАМЕТРЫ (только для расчёта масштабов)
# =========================================================
gamma = 2.211e5
alpha = 0.005
A_ex = 1.3e-11
Ms = 8e5
mu0 = 4 * np.pi * 1e-7
# =========================================================
# 2. МАСШТАБЫ
# l0 = sqrt(2A / (mu0*Ms^2)) → J = 1
# t0 = 1 / (gamma * Ms)
# =========================================================
l0 = np.sqrt(2.0 * A_ex / (mu0 * Ms**2)) # ≈ 5.69 нм
t0 = 1.0 / (gamma * Ms) # ≈ 5.65 пс
# =========================================================
# 3. БЕЗРАЗМЕРНЫЕ ВЕЛИЧИНЫ
# =========================================================
Nx, Ny = 64, 64
dx = 10e-9 / l0 # безразмерный шаг сетки ≈ 1.759
dt = 5e-13 / t0 # безразмерный шаг ≈ 0.0884
steps = 1200
save_every = 15
h_ext = 4.0e4 / Ms # безразмерное внешнее поле = 0.05
# Начальное условие: m = (1, 0, 0) везде
m = np.ones((Nx, Ny, 3)) * np.array([1.0, 0.0, 0.0])
# Возмущение в центре
cx, cy = Nx // 2, Ny // 2
angle = 0.52
radius_sq = 100
for ix in range(Nx):
for iy in range(Ny):
dist_sq = (ix - cx) ** 2 + (iy - cy) ** 2
if dist_sq < radius_sq:
m[ix, iy, 0] = np.cos(angle)
m[ix, iy, 1] = np.sin(angle)
# Нормировка один раз в начале
norm = np.sqrt(np.sum(m ** 2, axis=2, keepdims=True))
norm[norm == 0] = 1.0
m = m / norm
# =========================================================
# 4. БЕЗРАЗМЕРНОЕ ОБМЕННОЕ ПОЛЕ (J = 1)
# h_ex = nabla'^2 m (через roll, периодические границы)
# =========================================================
def compute_h(m_vec):
h = np.zeros_like(m_vec)
for i in range(3):
d2x = (np.roll(m_vec[:, :, i], -1, axis=0)
- 2 * m_vec[:, :, i]
+ np.roll(m_vec[:, :, i], 1, axis=0)) / dx**2
d2y = (np.roll(m_vec[:, :, i], -1, axis=1)
- 2 * m_vec[:, :, i]
+ np.roll(m_vec[:, :, i], 1, axis=1)) / dx**2
h[:, :, i] = d2x + d2y
return h
# =========================================================
# 5. БЕЗРАЗМЕРНОЕ УРАВНЕНИЕ ЛЛГ
# dm/dτ = -1/(1+α²) (m × h) - α/(1+α²) m × (m × h)
# =========================================================
def llg_rhs(m_vec):
h = compute_h(m_vec)
h[:, :, 0] += h_ext
cross_mh = np.cross(m_vec, h)
cross_m_mh = np.cross(m_vec, cross_mh)
return -(1.0 / (1 + alpha**2)) * cross_mh \
- (alpha / (1 + alpha**2)) * cross_m_mh
# =========================================================
# 6. ИНТЕГРИРОВАНИЕ (RK2 — метод средней точки, без нормировки)
# =========================================================
frames = []
print("Считаем кадры для анимации...")
for step in range(steps):
k1 = llg_rhs(m)
k2 = llg_rhs(m + 0.5 * dt * k1)
m += dt * k2
if np.any(~np.isfinite(m)):
print(f"⚠️ NaN на шаге {step} — уменьшите dt")
break
if step % save_every == 0 or step == steps - 1:
frames.append(m.copy())
print(f"Всего кадров: {len(frames)}")
# Проверка нормы
norm_check = np.sqrt(np.sum(m**2, axis=2))
print(f"|m|: min={norm_check.min():.6f}, max={norm_check.max():.6f} (без нормировки)")
# =========================================================
# 7. АНИМАЦИЯ (GIF через Pillow, без FuncAnimation.save)
# =========================================================
fig, ax = plt.subplots(figsize=(8, 8))
im = ax.imshow(m[:, :, 0], cmap='coolwarm', origin='lower', vmin=-1, vmax=1)
cbar = fig.colorbar(im, ax=ax)
cbar.set_label(r'$m_x$')
stride = 4
x_vec, y_vec = np.meshgrid(
np.arange(0, Nx, stride),
np.arange(0, Ny, stride),
indexing='ij'
)
quiver = ax.quiver(x_vec, y_vec,
np.zeros_like(x_vec), np.zeros_like(y_vec),
color='black', width=0.003, scale=30)
ax.set_title("2D спиновая волна: распространение возмущения")
ax.set_xlabel("x (ячейка)")
ax.set_ylabel("y (ячейка)")
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, fontsize=12,
bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
frames_pil = []
for idx in range(len(frames)):
m_frame = frames[idx]
t_ns = idx * save_every * dt * t0 * 1e9
im.set_data(m_frame[:, :, 0])
u = m_frame[::stride, ::stride, 0]
v = m_frame[::stride, ::stride, 1]
quiver.set_UVC(u, v)
time_text.set_text(f"t = {t_ns:.2f} нс")
buf = io.BytesIO()
fig.savefig(buf, format='png', dpi=100)
buf.seek(0)
frames_pil.append(Image.open(buf).convert('RGB'))
plt.close(fig)
gif_file = 'spin_wave_2d_nondim.gif'
if frames_pil:
frames_pil[0].save(
gif_file,
save_all=True,
append_images=frames_pil[1:],
duration=50,
loop=0
)
print(f"✅ GIF сохранён: {gif_file}")Сделал обезразмеривание в исходном 2D коде.Результат теперь принципиально другой:
Стало реалистичнее?
Проблема с пространственно разрывной намагниченностью решена?
Теперь стало гораздо лучше, но вопросы все равно остались. 1) В последних программах функции llg_rhs разные, хотя и не зависят от вычисления эффективного поля. Почему? 2) Во второй программе, шаг решетки в эффективных единицах 1.759, это достаточно малый? Сравните с шагом по времени. Как это сравнивается с шагом решетки в первой программе?
Расчет с доменной стенкой подозрительный. Разные участки должны вращаться с разной частотой (почему?). Это в конце концов должно приводить к расплыванию блоховской доменной стенки, потому что она не является собственным состоянием в магнетике с нулевой анизотропией. На картинке, однако, частота выглядит одинаковой и структура доменной стенки сохраняется.
Нужно помнить, что вера в код - это не религиозное чувство (очень хочется, чтобы код был правильным), а вынужденное состояние (нет ни единого сомнения в его правильности). В железобетонном научном исследовании действует презумпция виновности: любой код, любая формула считаются неверными, если не доказано обратное. Следовать ли этому принципу и если да, то насколько, каждый решает для себя сам в каждой конкретной ситуации. Но в целом старая поговорка "Нашел эффект - ищи дефект" дает хорошее общее правило.
1. Почему llg_rhs разные? Форма Гилберта ($1/(1+\alpha^2)$) vs форма Ландау–Лифшица (без неё) — разница пренебрежимо мала при $\alpha=0.005$. Но знак damping отличался: 1D-оригинал имел + (антидемпфинг при $\gamma > 0$), 2D-оригинал имел - (правильно). При обезразмеривании я поменял знак, не сказав об этом.
2. Шаг сетки. Без μ0: $dy' \approx 157$ — огромный, тривиально устойчивый, но физика неверна. С μ0: $dy' \approx 0.18$, и критерий устойчивости RK2 для LLG — не $dt \leq dy^2/2$ (это для диффузии), а $u = 4dt/dy^2 \lesssim 0.3$ (прецессия + слабый демпфинг). 2D-код: $u \approx 0.11$ — на грани, но устойчив.
3. Стенка не расплывалась из-за комбинации всех трёх багов. «Презумпция виновности» полностью оправдалась.
Вот исправленный код:
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io
# =========================================================
# 1. РАЗМЕРНЫЕ ПАРАМЕТРЫ
# =========================================================
Ms = 800e3; A_ex = 13e-12; alpha = 0.005; gamma = 2.211e5
mu0 = 4 * np.pi * 1e-7 # ← ВСЁ В SI, μ0 ОБЯЗАТЕЛЕН
# =========================================================
# 2. МАСШТАБЫ (с μ0!)
# l0 = sqrt(2A/(μ0·Ms²)) → J = 1
# t0 = 1/(γ·Ms)
# =========================================================
l0 = np.sqrt(2.0 * A_ex / (mu0 * Ms**2)) # ≈ 5.686 нм
t0 = 1.0 / (gamma * Ms) # ≈ 5.654 пс
# =========================================================
# 3. БЕЗРАЗМЕРНЫЕ ВЕЛИЧИНЫ
# =========================================================
J = 1.0
h_ext_x = 4e4 / Ms # ≈ 0.05
L = 200e-9 / l0 # ≈ 35.17
N = 101
y = np.linspace(-L/2, L/2, N)
dy = y[1] - y[0] # ≈ 0.352
# Устойчивость RK2 для LLG: u = 4·dt/dy² ≲ 0.3
dt = 0.003 # u ≈ 0.097 ✓
tau_max = 200.0 # ≈ 1.13 нс
n_steps = int(tau_max / dt)
save_every = max(1, n_steps // 50)
# --- Начальное условие: блоховская стенка ---
delta = 20e-9 / l0 # ≈ 3.52
theta = np.pi/2 * (1.0 + np.tanh(y / delta))
m = np.stack([np.sin(theta), np.cos(theta), np.zeros_like(y)], axis=0)
# =========================================================
# 4. ОБМЕННОЕ ПОЛЕ (J=1, безразмерное)
# =========================================================
def lap1d(f, dy):
return (np.roll(f, -1) - 2*f + np.roll(f, 1)) / dy**2
def compute_h(mv):
return np.stack([
lap1d(mv[0], dy) + h_ext_x,
lap1d(mv[1], dy),
lap1d(mv[2], dy)
], axis=0)
# =========================================================
# 5. УРАВНЕНИЕ ЛЛГ (Gilbert, ПРАВИЛЬНЫЙ ЗНАК damping = МИНУС)
# dm/dτ = -1/(1+α²)·(m×h) - α/(1+α²)·(m×(m×h))
# =========================================================
def llg_rhs(mv):
h = compute_h(mv)
c1 = np.cross(mv, h, axis=0)
c2 = np.cross(mv, c1, axis=0)
return -(1.0 / (1 + alpha**2)) * c1 \
- (alpha / (1 + alpha**2)) * c2
# =========================================================
# 6. ИНТЕГРИРОВАНИЕ (RK2, БЕЗ нормировки)
# =========================================================
taus, my_p, mx_p, mz_p = [], [], [], []
mc = m.copy()
for step in range(n_steps):
k1 = llg_rhs(mc)
k2 = llg_rhs(mc + 0.5 * dt * k1)
mc += dt * k2
if step % save_every == 0 or step == n_steps - 1:
taus.append(step * dt)
my_p.append(mc[1].copy())
mx_p.append(mc[0].copy())
mz_p.append(mc[2].copy())
taus = np.array(taus)
my_p = np.array(my_p)
mx_p = np.array(mx_p)
mz_p = np.array(mz_p)
times_real = taus * t0
y_nm = y * l0 * 1e9
norm = np.sqrt(mc[0]**2 + mc[1]**2 + mc[2]**2)
print(f"l0={l0*1e9:.3f} нм, t0={t0*1e12:.3f} пс, J={J}")
print(f"dy={dy:.4f}, dt={dt:.4f}, u={4*dt/dy**2:.4f}")
print(f"Шагов: {n_steps}, кадров: {len(taus)}")
print(f"|m|: min={norm.min():.6f}, max={norm.max():.6f} (без нормировки)")
print(f"Время до {times_real[-1]*1e9:.2f} нс")
# =========================================================
# 7. ГРАФИК: 3 панели (m_y, m_z, |dm/dy|)
# =========================================================
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
indices = [0, len(taus)//4, len(taus)//2, 3*len(taus)//4, -1]
colors = ['#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd']
for i, idx in enumerate(indices):
t_ns = times_real[idx] * 1e9
axes[0].plot(y_nm, my_p[idx], color=colors[i], lw=2, label=f't={t_ns:.2f} нс')
axes[1].plot(y_nm, mz_p[idx], color=colors[i], lw=2, label=f't={t_ns:.2f} нс')
dmx = np.gradient(mx_p[idx], y_nm)
dmy = np.gradient(my_p[idx], y_nm)
dmz = np.gradient(mz_p[idx], y_nm)
axes[2].plot(y_nm, np.sqrt(dmx**2+dmy**2+dmz**2), color=colors[i], lw=2, label=f't={t_ns:.2f} нс')
axes[0].set_title(r'$m_y$'); axes[0].set_ylabel(r'$m_y$')
axes[1].set_title(r'$m_z$ (3D прецессия!)'); axes[1].set_ylabel(r'$m_z$')
axes[2].set_title(r'$|d\mathbf{m}/dy|$ (расплывание)'); axes[2].set_ylabel('нм$^{-1}$')
for ax in axes:
ax.set_xlabel('y (нм)'); ax.legend(fontsize=9); ax.grid(True, ls=':', alpha=0.3)
plt.tight_layout()
plt.savefig('wall_analysis.png', dpi=150)
plt.close()
# =========================================================
# 8. GIF
# =========================================================
fig2, ax = plt.subplots(figsize=(10, 6))
line, = ax.plot([], [], lw=2, color='#1f77b4')
ax.set_xlim(y_nm[0], y_nm[-1]); ax.set_ylim(-1.1, 1.1)
ax.axhline(0, color='k', ls='--', lw=1, alpha=0.4)
ax.set_xlabel('y (нм)'); ax.set_ylabel(r'$m_y$')
title = ax.set_title(''); ax.grid(True, ls=':', alpha=0.3)
frames = []
for i in range(len(taus)):
line.set_data(y_nm, my_p[i])
title.set_text(f't = {times_real[i]*1e9:.2f} нс')
buf = io.BytesIO()
fig2.savefig(buf, format='png', dpi=80)
buf.seek(0)
frames.append(Image.open(buf).convert('RGB'))
plt.close(fig2)
frames[0].save('wall_spreading.gif', save_all=True,
append_images=frames[1:], duration=50, loop=0)
print("✅ wall_analysis.png, wall_spreading.gif")
В результате получим:
Теперь всё верно? Дан ответ на все вопросы? Проблема решена? Полученная модель соответствует реальности?
Я бы и 
во временной масштаб отправил, но в интересных случаях слабого затухания это, действительно, непринципиально, просто формулы микроскопически лучше выглядят.
Кошмар с доменной стенкой происходит от наложенных периодических условий. Я в этом отношении ненастоящий теоретический физик: периодические граничные условия мне не нравятся. Здесь можно проговорить мантру про самосопряженные расширения оператора импульса, но сейчас это неважно. В данном случае, это приводит к эффективной разрывной намагниченности и мгновенному скачку намагниченности за один шаг по времени. Как только такое происходит на все остальное можно махнуть рукой. В частности поэтому в своем самом первом коментарии я упомянул физичность начальных условий.
А на вопрос соответствует ли полученная модель (даже с полностью выправленными деталями) реальности дам тот же ответ, что в свое время получил сам "я не знаю". Реальность - штука сложная. Там будут размагничивающие факторы, свои явления на границах и не будет двумерных систем.
А как можно исправить исходную статью? Полностью переделывать или писать новую длинную на ту же тему я не собираюсь.
Исходную статью никак не исправить. Она содержит характерную фразу
Наблюдаемая картина согласуется с теоретическими представлениями о спиновых волнах в ферромагнетиках.
И это при том, что на всей этой странице не появилось ни одной картинки, которая бы показывала что-то большее, чем артефакты неправильно организованного численного счета. Это делает процитированную фразу довольно забавной и показывает ценность таких высказываний.
Статьи, вообще, имеет смысл писать только о том, что важно для самого себя: удалось в чем-то разобраться, получилось увидеть что-то новое в, казалось бы, обычной ситуации и т.д.
А ведь вы тоже защитили кандидатскую по теме спиновой динамики?
Спиновые волны были одной из составляющих.
В ходе своей научной работы вы делали что-то подобное? То есть проводили моделирование уравнение ЛЛГ, писали код, получали подобные картинки, гифки?
Все перечисленное плюс много аналитического счета. Для меня именно аналитика является главным источником знаний.
Можете поделиться своими картинками? А экспериментировали ли вы с реальными ферромагнетиками (пермаллоем, и т. д.)?
Главным образом картинки у меня в публикациях, которые можно в Google Scholar посмотреть. В рабочем порядке картинок генерируется порядка на два-три больше и отношение к ним как к бросовому материалу. Первопринципным моделированием ЛЛГ не занимался больше четверти века, да и тогда немного. Тех картинок в принципе не разыскать.
С ферромагнетиками не экспериментировал.
Как это происходит в реальности, наверное, знает Кирилл Циберкин. Он по этой теме защитил докторскую диссертацию и последние 5 лет работал над этой задачей.
Я добавил ещё оценку устойчивости при помощи показателя Ляпунова, чтобы обосновать это всё:
Устойчивость динамики намагниченности оценивалась по локальному показателю Ляпунова, рассчитанному численно путём сравнения двух близких траекторий в фазовом пространстве. При малых значениях внешнего поля показатель λ<0, что соответствует устойчивым регулярным колебаниям (режим ускорения). В промежуточной области наблюдается переход к λ>0, указывающий на возникновение хаотической динамики и чувствительность к начальным условиям. При высоких полях показатель снова становится отрицательным (λ<0), что свидетельствует о восстановлении устойчивости и переходе к режиму вынужденной синхронизации.
При этом показатель Ляпунова на всём диапазоне колебался около нуля, но на подавляющем большинстве интервалов оставался отрицательным.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Спиновая динамика в анимациях: моделируем волны намагниченности в ферромагнетиках с помощью python