Численные методы решения дифференциальных уравнений

@Tatiana_L · Регистрация: 25.08.2013

Студворк — интернет-сервис помощи студентам

Дамы и господа!
Имеется дифференциальное уравнения типа d^2y/dx^2=a*exp(-b*x)-c*(dy/dx)^2/y-n*(dy/dx)/x-m*y^3,
где a,b,c,n,m -- числа, константы.
Нужно решить численными методами это уравнения с условием того, что имеются начальные условия.
Однако, решая это уравнение классическим методом Рунге-Кутта 4-порядка, возникает сингулярность при переходе y через 0. Вероятно потому, что при делении на очень маленькое число, одно из слагаемых начинает стремительно возрастать.
В связи с этим у меня вопрос: возможно кто-то сталкивался уже с таким и знает как обойти сингулярность и решить данное уравнение численным методом или предложит численный метод для решения данного уравнения?

@Том Ардер · 25.08.2013, 21:15

Комментарий модератора

Tatiana_L, Правила, 4.7. Как можно более полно описывайте суть проблемы или вопроса, что было сделано для ее решения и какие результаты получены.
Редактор формул

@Tatiana_L · 26.08.2013, 01:09 **[ТС]**

уравнение $https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{\operatorname{d}^{2}y}{\operatorname{d}{x}^{2}}=12*\pi *exp{-2*x}-\frac{1}{2*y}*(\frac{\operatorname{d}y}{\operatorname{d}x})^2-\frac{2}{x}*(\frac{\operatorname{d}y}{\operatorname{d}x})-4*y^3$ .
Начальные условия для решения данного уравнения $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(0.35)=1.8145, \frac{\operatorname{d}y}{\operatorname{d}x}(0.35)=-1.4804$
Последняя точка, которая была посчитана $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(3.84)=0.00573$ .
Результаты расчетов находятся в excel файле "уравнение".

@cmath · 26.08.2013, 01:39

Tatiana_L, на каком отрезке вы интегрируете это уравнение?

@Tatiana_L · 26.08.2013, 01:41 **[ТС]**

Мне нужно как можно ближе к 0 и с правой стороны как мин до 10, но было бы лучше, если бы и дальше.
т е задача минимум (0,10].

@Том Ардер · 27.08.2013, 01:05

Предположим, что в нуле решение ведёт себя степенным образом:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?y(x)\sim cx^\alpha, \; x\rightarrow 0$
Подстановка в исходное уравнение даёт два варианта: регулярное решение при $https://www.cyberforum.ru/cgi-bin/latex.cgi?\alpha$ = 2
и сингулярное при -1 < $https://www.cyberforum.ru/cgi-bin/latex.cgi?\alpha$ < 0.
Ввиду нелинейности уравнения коэффициент c зависит от $https://www.cyberforum.ru/cgi-bin/latex.cgi?\alpha$ .
Выбираем регулярное решение, делаем подстановку $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(x)=cx^2z(x), \; z(0)=1$
и решаем уравнение для новой функции z(x).

240Volt · 27.08.2013, 02:38

Tatiana_L! Давайте посмотрим на правую часть вашего уравнения.
В окрестности особой точки (где-то около полученного вами 3,84) в правой части уравнения главную роль играет второе слагаемое.
Поэтому решение вашего уравнения в узкой окрестности этой точки имеет степенную особенность вида
$https://www.cyberforum.ru/cgi-bin/latex.cgi?y \approx A{\left| {b - x} \right|^{2/3}}$ .
Можно было бы попробовать выделить эту особенность, записав решение в виде
$https://www.cyberforum.ru/cgi-bin/latex.cgi?y = {\left| {b - x} \right|^{2/3}}u(x)$ .
Тогда уравнение для регулярной части u(x) уже не имело бы особых точек, кроме нуля, и с вашими начальными условиями (не в нуле) вы это уравнение решили бы.
Проблема, однако, в том, что особая точка b неизвестна, а с неправильным b уравнение для u(x) все равно не удастся решить численно.
Можно, конечно, попробовать взять от безысходности $https://www.cyberforum.ru/cgi-bin/latex.cgi?b = 3,84$ , но ничего хорошего из этого, скорее всего, не выйдет.
Я бы попробовал решать уравнение для u(x) с разными значениями b, пока не найдется такое значение b, для которого сингулярность в численном решении не возникает, либо находится далеко за пределами нужной вам области значений x.

@Том Ардер · 27.08.2013, 14:24

240Volt, навёл на идею. Удалось подобрать замену: $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(x)=(u(x) )^\nu$
При $https://www.cyberforum.ru/cgi-bin/latex.cgi?^\nu=2/3$ сингулярное слагаемое уходит из уравнения. Однако решение заметно отличается от исходного. Причины пока непонятны. Файл Mathcad (вер. 14) прикрепляю.

@Том Ардер · 27.08.2013, 20:14

Исправил ошибку в коэффициенте, решения совпали.

@Tatiana_L · 28.08.2013, 09:13 **[ТС]**

Спасибо!
это очень помогло!
теперь осталось понять какой именно алгоритм решения использует Mathcad (дело в том, что я использую Maple, а он по умолчанию решает rkf_45, и в данном случае преобразование не работает там, но в Matcad все ок, я перепроверяла.)
Хочу еще это реализовать в excel.

@Том Ардер · 28.08.2013, 11:50

Tatiana_L, Mathcad по умолчанию использует метод Адамса-Башфорта. В контекстном меню функции Odesolve есть выбор других методов - Рунге-Кутта (4-го порядка) с фиксированным или адаптивным шагом.

Сообщение от Tatiana_L

в данном случае преобразование не работает там

Непонятно. Что может не нравиться Maple в преобразованном уравнении:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{d^2u}{dx^2}=18\pi \exp (-2x){u}^{1/3}-\frac{2}{x}\cdot \frac{du}{dx}-6{u}^{7/3}$

@Tatiana_L · 28.08.2013, 17:08 **[ТС]**

Да, у меня получилось точно такое уравнение.
Точки до "чего-то" лежат точно на точках, полученных при решении оригинального уравнения.Но дальше -- снова не хочет, причем как в excel, так и, как я уже и говорила, в Maple.

@Tatiana_L · 28.08.2013, 18:08 **[ТС]**

прикрепляю maple-файл с rkf_45 -- сингулярность остается(это в подтверждение своих слов)).

@Tatiana_L · 28.08.2013, 18:46 **[ТС]**

Сообщение от Том Ардер

Tatiana_L,
Непонятно. Что может не нравиться Maple в преобразованном уравнении:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{d^2u}{dx^2}=18\pi \exp (-2x){u}^{1/3}-\frac{2}{x}\cdot \frac{du}{dx}-6{u}^{7/3}$

это все в подтверждения вышесказанного.(скриншот на случай, если нет maple.

@Том Ардер · 28.08.2013, 20:35

Похоже, что причина в способах вычисления кубического корня и степени 1/3 - в математических пакетах они разные. Когда в Mathcad в уравнении заменил кубический корень на степень 1/3 - тоже перестал считать, и даже не доходя до 3.84, а немного раньше, при 3.68. Предполагаю аналогичную фичу (или баг) и в Maple (пакета нет, попробую развернуть Maple Player).

@Том Ардер · 28.08.2013, 23:00

Применил Wolfram Mathematica. Там тоже с кубическим корнем и степенью 1/3 своеобразие: после х = 3.84 в решении появляется небольшая мнимая часть. А при х = 9.28 происходит численная катастрофа. Графики реальной и мнимой части решения прилагаю.

@Tatiana_L · 28.08.2013, 23:05 **[ТС]**

Спасибо!
получилось исправить данную проблему в Maple(файл и скриншот прикрепляю), однако пока не удается в excel, так как там нет функции кубического корня(хотя есть sqrt).

@Tatiana_L · 28.08.2013, 23:12 **[ТС]**

Чтото похожее на второй график у меня получилось, когда я пыталась както упростить уравнение по второй части -- т е взять производную с числителя и знаменателя(часть, где первая производная в квадрате).

@Том Ардер · 29.08.2013, 11:33

Сообщение от Tatiana_L

пока не удается в excel, так как там нет функции кубического корня

Есть функция POWER(x;a), вычисляющая любые степени х^а. С отрицательными основаниями и степенью 1/3 работает правильно.

В Mathematica сделал правильную обработку кубического корня, теперь считает везде, и результаты совпадают с другими решениями.

@Tatiana_L · 30.08.2013, 04:48 **[ТС]**

Да, так я тоже пробовала ранее, но это не помогло.
думаю, что дело в самом алгоритме, потому что в excel я делала просто Рунге-Кутта 4.

Новые блоги и статьи Все статьи Все блоги /
Модульная разработка через nuget packages DevAlt 07.03.2026 Сложившийся в . Net-среде способ разработки чаще всего предполагает монорепозиторий в котором находятся все исходники. При создании нового решения, мы просто добавляем нужные проекты и имеем. . .	Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .
SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip Сканируйте QR-код на мобильном и вы увидите, что появится джойстик для управления главным героем. . . .	Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .	SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .

Численные методы решения дифференциальных уравнений

Решение

Решение

@Tatiana_L 0 / 0 / 0 Регистрация: 25.08.2013 Сообщений: 10

	Численные методы решения дифференциальных уравнений 25.08.2013, 19:54. Показов 3936. Ответов 20 Метки нет (Все метки) Дамы и господа! Имеется дифференциальное уравнения типа d^2y/dx^2=aexp(-bx)-c(dy/dx)^2/y-n(dy/dx)/x-m*y^3, где a,b,c,n,m -- числа, константы. Нужно решить численными методами это уравнения с условием того, что имеются начальные условия. Однако, решая это уравнение классическим методом Рунге-Кутта 4-порядка, возникает сингулярность при переходе y через 0. Вероятно потому, что при делении на очень маленькое число, одно из слагаемых начинает стремительно возрастать. В связи с этим у меня вопрос: возможно кто-то сталкивался уже с таким и знает как обойти сингулярность и решить данное уравнение численным методом или предложит численный метод для решения данного уравнения? 0

@cmath 2525 / 1751 / 152 Регистрация: 11.08.2012 Сообщений: 3,349
	26.08.2013, 01:39
	Tatiana_L, на каком отрезке вы интегрируете это уравнение? 1

@Tatiana_L 0 / 0 / 0 Регистрация: 25.08.2013 Сообщений: 10
	26.08.2013, 01:41 [ТС]
	Мне нужно как можно ближе к 0 и с правой стороны как мин до 10, но было бы лучше, если бы и дальше. т е задача минимум (0,10]. 0

@Том Ардер $Эксперт по математике/физике$ 4218 / 3418 / 396 Регистрация: 15.06.2009 Сообщений: 5,818
	27.08.2013, 01:05
	Предположим, что в нуле решение ведёт себя степенным образом: $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(x)\sim cx^\alpha, \; x\rightarrow 0$ Подстановка в исходное уравнение даёт два варианта: регулярное решение при $https://www.cyberforum.ru/cgi-bin/latex.cgi?\alpha$ = 2 и сингулярное при -1 < $https://www.cyberforum.ru/cgi-bin/latex.cgi?\alpha$ < 0. Ввиду нелинейности уравнения коэффициент c зависит от $https://www.cyberforum.ru/cgi-bin/latex.cgi?\alpha$ . Выбираем регулярное решение, делаем подстановку $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(x)=cx^2z(x), \; z(0)=1$ и решаем уравнение для новой функции z(x). 2

240Volt 4445 / 2449 / 227 Регистрация: 20.08.2011 Сообщений: 3,108
	27.08.2013, 02:38
	Сообщение было отмечено как решение Решение Tatiana_L! Давайте посмотрим на правую часть вашего уравнения. В окрестности особой точки (где-то около полученного вами 3,84) в правой части уравнения главную роль играет второе слагаемое. Поэтому решение вашего уравнения в узкой окрестности этой точки имеет степенную особенность вида $https://www.cyberforum.ru/cgi-bin/latex.cgi?y \approx A{\left\| {b - x} \right\|^{2/3}}$ . Можно было бы попробовать выделить эту особенность, записав решение в виде $https://www.cyberforum.ru/cgi-bin/latex.cgi?y = {\left\| {b - x} \right\|^{2/3}}u(x)$ . Тогда уравнение для регулярной части u(x) уже не имело бы особых точек, кроме нуля, и с вашими начальными условиями (не в нуле) вы это уравнение решили бы. Проблема, однако, в том, что особая точка b неизвестна, а с неправильным b уравнение для u(x) все равно не удастся решить численно. Можно, конечно, попробовать взять от безысходности $https://www.cyberforum.ru/cgi-bin/latex.cgi?b = 3,84$ , но ничего хорошего из этого, скорее всего, не выйдет. Я бы попробовал решать уравнение для u(x) с разными значениями b, пока не найдется такое значение b, для которого сингулярность в численном решении не возникает, либо находится далеко за пределами нужной вам области значений x. 3

@Том Ардер $Эксперт по математике/физике$ 4218 / 3418 / 396 Регистрация: 15.06.2009 Сообщений: 5,818
	27.08.2013, 20:14
	Исправил ошибку в коэффициенте, решения совпали. Миниатюры 2

@Tatiana_L 0 / 0 / 0 Регистрация: 25.08.2013 Сообщений: 10
	28.08.2013, 09:13 [ТС]
	Спасибо! это очень помогло! теперь осталось понять какой именно алгоритм решения использует Mathcad (дело в том, что я использую Maple, а он по умолчанию решает rkf_45, и в данном случае преобразование не работает там, но в Matcad все ок, я перепроверяла.) Хочу еще это реализовать в excel. 0

@Tatiana_L 0 / 0 / 0 Регистрация: 25.08.2013 Сообщений: 10
	28.08.2013, 17:08 [ТС]
	Да, у меня получилось точно такое уравнение. Точки до "чего-то" лежат точно на точках, полученных при решении оригинального уравнения.Но дальше -- снова не хочет, причем как в excel, так и, как я уже и говорила, в Maple. 0

@Том Ардер $Эксперт по математике/физике$ 4218 / 3418 / 396 Регистрация: 15.06.2009 Сообщений: 5,818
	28.08.2013, 20:35
	Похоже, что причина в способах вычисления кубического корня и степени 1/3 - в математических пакетах они разные. Когда в Mathcad в уравнении заменил кубический корень на степень 1/3 - тоже перестал считать, и даже не доходя до 3.84, а немного раньше, при 3.68. Предполагаю аналогичную фичу (или баг) и в Maple (пакета нет, попробую развернуть Maple Player). 2

@Том Ардер $Эксперт по математике/физике$ 4218 / 3418 / 396 Регистрация: 15.06.2009 Сообщений: 5,818
	28.08.2013, 23:00
	Применил Wolfram Mathematica. Там тоже с кубическим корнем и степенью 1/3 своеобразие: после х = 3.84 в решении появляется небольшая мнимая часть. А при х = 9.28 происходит численная катастрофа. Графики реальной и мнимой части решения прилагаю. Миниатюры 1

@Tatiana_L 0 / 0 / 0 Регистрация: 25.08.2013 Сообщений: 10
	28.08.2013, 23:12 [ТС]
	Чтото похожее на второй график у меня получилось, когда я пыталась както упростить уравнение по второй части -- т е взять производную с числителя и знаменателя(часть, где первая производная в квадрате). 0

@Tatiana_L 0 / 0 / 0 Регистрация: 25.08.2013 Сообщений: 10
	30.08.2013, 04:48 [ТС]
	Да, так я тоже пробовала ранее, но это не помогло. думаю, что дело в самом алгоритме, потому что в excel я делала просто Рунге-Кутта 4. 0