Не решается диффур через ode 45 в матлаб 2017b

[optimden]

Здравствуйте. Возникла проблема с решением диффура в матлабе через ode45.
На следующий код матлаб 2017b долго не реагирует, не выводит массив n1(r), не строит график:

Используется синтаксис Matlab M

lambda = 3.8709*10^-9;
diap = [0.00001 0.47];
n00=[10^17 0];
[r1,n1]=ode45(@(r1,n1)[n1(2);-(1/r1).*n1(2)-(1/(lambda).^2).*n1(1)],diap,n00);
disp(n1(:,1));
plot(r1(:,1),n1(:,1));

Ошибок при этом никаких не пишет.
Как можно решить данную проблему?
Заранее признателен.

[Pphantom]

Если я правильно восстановил исходную задачу (я не любитель MATLAB'а, так что мог и ошибиться), вы пытаетесь решить уравнение вида

$$ \ddot x + \frac{\dot x}{t} + \frac{x}{\lambda} = 0 $$

для $x(t)$ на интервале $t \in [10^{-5}, 0.47]$ с начальными условиями $x(t_0) = 10^{17}$ и $\dot x(t_0) = 0$ (где $t_0 = 10^{-5}$). Параметр $\lambda = 3.8709 \cdot 10^{-9}$.

Отсюда вопрос: вы про вычислительную математику что-нибудь слышали?

В общем-то то, что вы пытаетесь сделать с MATLAB'ом - это такое своеобразное изнасилование. Понадеяться на то, что явные методы Рунге-Кутты 4 или 5 порядка (а ode45 реализует именно их) справятся с существенно жесткой системой, конечно, можно, но практически интегратор с адаптивным шагом будет делать этот шаг все меньше и меньше, и результат вы получите (если получите) очень и очень нескоро.

Так что это не проблема MATLAB'а, он-то все делает правильно, просто вы хотите от него слишком многого. Если эта задача имеет какую-то предшествующую историю, излагайте: возможно, ее можно переформулировать иначе. Если вам действительно нужно решить это уравнение именно в таком виде - нужно подбирать метод решения с учетом ее особенностей.

[optimden]

Попробовал сейчас решить символьным методом:

Используется синтаксис Matlab M

syms n(r) r;
equ = diff(n,r,2)+(1/r)*diff(n,r)==n/(lambda)^2;
Dn = diff(n,r);
cond = ([n(0.00001)==10^17 Dn(0.00001)==1]);
nr(r) = dsolve(equ,cond);
r = 0.00001:0.01:0.47;
plot(nr(r));

В таком случае решение матлаб находит в виде комбинации функций Бесселя, но график тоже почему-то не строит.

-- 19.05.2019, 15:10 --

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

$$\ddot n+\frac{\dot n}{r}+\frac{n}{\lambda}=0$$

Решением является функция Бесселя нулевого порядка.

$$\lambda=\frac{D_a}{v_i}$$

где $D_a$ - коэффициент амбиполярной диффузии, $v_i$ - частота ионизации.

[GAA]

При заданных значениях параметров на начальном участке решение возможно быстро возрастает.

Вложение:

n.PNG [ 7.63 Кб | Просмотров: 1886 ]

Попробуйте уменьшить промежуток.

-- Sun 19.05.2019 14:50:52 --

Так можно на скорую руку график построить

Используется синтаксис Matlab M

r = 0.00001:0.0000000001:0.0000101;
Y= double(vpa(nr(r)))
plot(r, Y);

[Pphantom]

optimden в сообщении #1394011 писал(а):

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

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

[optimden]

Нормировал концентрацию, сделал $\lambda=1$:

Используется синтаксис Matlab M

syms n(r) r;
equ = diff(n,r,2)+(1/r)*diff(n,r)==n;
Dn = diff(n,r);
cond = ([n(0.000001)==1 Dn(0.000001)==0]);
nr(r) = dsolve(equ,cond);
r = 0.000001:0.0001:0.0065;
disp(nr(r));
plot(r,nr(r));

Также изменил диапазон $r$. Всё равно матлаб не выдаёт ни графика, ни ошибок.

[GAA]

Используется синтаксис Matlab M

syms n(r) r;
equ = diff(n,r,2)+(1/r)*diff(n,r)==n;
Dn = diff(n,r);
cond = ([n(0.000001)==1 Dn(0.000001)==0]);
nr(r) = dsolve(equ,cond);
r = 0.000001:0.0001:0.0065;
Y= vpa(nr(r));
plot(r, Y);

Вложение:

Комментарий к файлу: matlab 2013b

n.PNG [ 9.62 Кб | Просмотров: 1855 ]

-- Sun 19.05.2019 17:11:05 --

r = 0.000001:0.0001:0.0065;
Теперь промежуток можно увеличить. Однако замену независимого переменного Вы, на мой взгляд, выполнили неправильно или в исходном тексте не все изменили.

[optimden]

А что означает

Используется синтаксис Matlab M

Y = vpa(nr(r));

?
То, что график получается, уже неплохо, но в итоге должна получаться функция Бесселя нулевого рода, что и видно при выводе символьного решения через $disp\left(nr\left(r\right)\right);$
А график тем не менее совсем не похож на график функции Бесселя.

-- 19.05.2019, 18:14 --

Вы про замену $lambda$ на 1?

[GAA]

optimden в сообщении #1394049 писал(а):

А что означает

Уже писал: post1387110.html#p1387110.

У меня R2013b возвращает линейную комбинацию функций Бесселя первого и второго рода (besseli(0, r), besselk(0, r) ):

nr(r) =
 (besseli(1, 1/1000000)*besselk(0, r))/(besseli(0, 1/1000000)*besselk(1, 1/1000000) + besseli(1, 1/1000000)*besselk(0, 1/1000000)) + (besselk(1, 1/1000000)*besseli(0, r))/(besseli(0, 1/1000000)*besselk(1, 1/1000000) + besseli(1, 1/1000000)*besselk(0, 1/1000000))

(Нулевого рода/типа я функций Бесселя не встречал.)

-- Sun 19.05.2019 17:34:17 --

optimden в сообщении #1394049 писал(а):

Вы про замену $\lambda$ на 1?

Надо же не просто положить $\lambda=1$. Надо выполнить замену независимого переменного. В начальном сообщении $\lambda$ была в квадрате, тогда замена: $r=\lambda \rho$.

-- Sun 19.05.2019 17:37:28 --

В более компактном виде решение (с приближёнными значениями коэффициентов):

>> digits(12)
>> vpa(nr)
ans(r) = 0.999999999993*besseli(0.0, r) + 5.0e-13*besselk(0.0, r)

-- Sun 19.05.2019 18:15:10 --

Вместо plot можно попробовать воспользоваться устаревшей ezplot:

Используется синтаксис Matlab M

syms n(r) r;
equ = diff(n,r,2)+(1/r)*diff(n,r)==n;
Dn = diff(n,r);
cond = ([n(0.000001)==1 Dn(0.000001)==0]);
nr(r) = dsolve(equ,cond)
ezplot(nr, [0.000001 0.0065])
title('nr')

R2013b рисует.

upd. Графики модифицированных функций Бесселя можно посмотреть в статье Википедии. (Для нулевого порядка на глаз соответствует.)

[optimden]

Большое спасибо всем ответившим. Задача в итоге решилась через ode45.