2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 13:08 
Здравствуйте. Возникла проблема с решением диффура в матлабе через 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));

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

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 14:54 
Если я правильно восстановил исходную задачу (я не любитель 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}$.

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

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

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

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 15:01 
Попробовал сейчас решить символьным методом:
Используется синтаксис 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$ - частота ионизации.

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 15:39 
При заданных значениях параметров на начальном участке решение возможно быстро возрастает.
Вложение:
n.PNG

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

-- Sun 19.05.2019 14:50:52 --

Так можно на скорую руку график построить
Используется синтаксис Matlab M
r = 0.00001:0.0000000001:0.0000101;
Y= double(vpa(nr(r)))
plot(r, Y);


У вас нет доступа для просмотра вложений в этом сообщении.

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 16:08 
optimden в сообщении #1394011 писал(а):
Да, вы верно интерпретировали задачу. Надо решить уравнение распределения концентрации электронов в плазме положительного столба от радиальной координаты трубки, в которой находится плазма:
Ну тогда первое действие, которое нужно было сделать - обезразмерить задачу, чтобы не считать результаты с гигантскими порядками. Ну и в силу физики задачи имело смысл считать не концентрацию, а ее логарифм.

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 17:29 
Нормировал концентрацию, сделал $\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$. Всё равно матлаб не выдаёт ни графика, ни ошибок.

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 18:06 
Используется синтаксис 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);
Вложение:
n.PNG


-- Sun 19.05.2019 17:11:05 --

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


У вас нет доступа для просмотра вложений в этом сообщении.

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 18:12 
А что означает
Используется синтаксис Matlab M
Y = vpa(nr(r));
?
То, что график получается, уже неплохо, но в итоге должна получаться функция Бесселя нулевого рода, что и видно при выводе символьного решения через $disp\left(nr\left(r\right)\right);$
А график тем не менее совсем не похож на график функции Бесселя.

-- 19.05.2019, 18:14 --

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

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 18:27 
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. Графики модифицированных функций Бесселя можно посмотреть в статье Википедии. (Для нулевого порядка на глаз соответствует.)

 
 
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение21.05.2019, 00:10 
Большое спасибо всем ответившим. Задача в итоге решилась через ode45.

 
 
 [ Сообщений: 10 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group