2014 dxdy logo

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

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




Начать новую тему Ответить на тему
 
 Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 13:08 


21/02/19
108
Здравствуйте. Возникла проблема с решением диффура в матлабе через 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 
Заслуженный участник


09/05/12
25179
Если я правильно восстановил исходную задачу (я не любитель 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 


21/02/19
108
Попробовал сейчас решить символьным методом:
Используется синтаксис 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 
Заслуженный участник


12/07/07
4438
При заданных значениях параметров на начальном участке решение возможно быстро возрастает.
Вложение:
n.PNG
n.PNG [ 7.63 Кб | Просмотров: 1646 ]

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

-- 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 
Заслуженный участник


09/05/12
25179
optimden в сообщении #1394011 писал(а):
Да, вы верно интерпретировали задачу. Надо решить уравнение распределения концентрации электронов в плазме положительного столба от радиальной координаты трубки, в которой находится плазма:
Ну тогда первое действие, которое нужно было сделать - обезразмерить задачу, чтобы не считать результаты с гигантскими порядками. Ну и в силу физики задачи имело смысл считать не концентрацию, а ее логарифм.

 Профиль  
                  
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 17:29 


21/02/19
108
Нормировал концентрацию, сделал $\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 
Заслуженный участник


12/07/07
4438
Используется синтаксис 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
n.PNG [ 9.62 Кб | Просмотров: 1615 ]


-- Sun 19.05.2019 17:11:05 --

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

 Профиль  
                  
 
 Re: Не решается диффур через ode 45 в матлаб 2017b
Сообщение19.05.2019, 18:12 


21/02/19
108
А что означает
Используется синтаксис 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 
Заслуженный участник


12/07/07
4438
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 


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 10 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group