Розв’язування
Ця задача, в
порівнянні з попередньою, має дві
принципові особливості. По перше, шукана
функція залежить тільки від двох змінних
– часу t
і відстані до центру мембрани
.
Тому для розв’язування задачі необхідно
перейти до полярної системи координат
(в цій системі координат незалежними
змінними є відстань до центра
і кут повороту
,
але від нього шукана функція не залежить.
По друге, всі початкові умови задачі
нульові, а граничні умови не стаціонарні
(залежать від часу). Тому необхідно
змінити алгоритм пошуку розв’язування.
Запишемо оператор
Лапласа в полярній системі координат
(в декартовій він має вигляд
)
врахувавши, що функція на яку він діє
залежить тільки від
і не залежить від
.
В цьому випадку корисна процедура
dchange().
> PDEtools[dchange]({x=rho*cos(phi),y=rho*sin(phi)},
diff(z(sqrt(x^2+y^2)),x$2)+diff(z(sqrt(x^2+y^2)), y$2),{phi,rho});
Одержаний вираз необхідно спростити.
> simplify(%,symbolic);
Тепер можна записати рівняння.
> Eqn:=diff(u(rho,t),t$2)=a^2*(1/rho)*diff(rho*diff(u(rho, t),rho),rho);
Крім того, задамо функціональну залежність, відповідно до якої відбувається рух країв мембрани.
> f:=t->A*sin(omega*t);
Спочатку знайдемо розв’язок, який буде задовольняти нестаціонарні граничні умови. Розв’язок будемо шукати методом розділення змінних.
> pdsolve(Eqn,HINT=F(rho)*f(t));
Функція
має задовольняти такі умови:
,
щоб забезпечити потрібний закон руху
границі мембрани, і, крім того, вона
повинна бути обмеженою нулі (F(0)<>infinity).
> dsolve({diff(F(rho),`$`(rho,2))=-F(rho)*omega^2*rho+
+a^2*diff(F(rho),rho))/a^2/rho,F(L)=1,F(0)<>infinity}, F(rho));
> F:=unapply(rhs(%),rho);
Тепер в початковому
рівнянні перейдемо до функції
,
яка пов’язана з початковою функцією
співвідношенням
,
і одержимо наступне.
> Eqn2:=subs(u(rho,t)=v(rho,t)+F(rho)*f(t),Eqn);
Дане рівняння спрощуємо.
> Eqn2:=simplify(Eqn2);
> Eqn2:=lhs(Eqn2)-rhs(Eqn2)=0;
> Eqn2:=simplify(Eqn2);
Розв’язуємо рівняння методом розділення змінних.
> pdsolve(Eqn,HINT=R(rho)*T(t));
Оскільки функція
обмежена
в нулі, то одержимо.
> dsolve({diff(R(rho),`$`(rho,2))=R(rho)*(-lambda^2)-diff(R(rho),rho)/rho,R(0)<>infinity},R(rho));
Параметр
має бути таким, щоб автоматично
виконувались граничні умови для функції
.
Ці умови нульові. Звідси одержимо задачу
на власні значення.
> solve(BesselJ(0,lambda*L)=0,lambda);
RootOf(BesselJ(0,_ZL))
Поданий в області
вводу результат означає, що власні
значення даної задачі одержуються
шляхом ділення на L
нульової функції Бесселя (маються на
увазі розв’язки рівняння
).
В Maple є процедура BesselJZeros(), яка дозволяє обчислити значення аргументу, при яких функція Бесселя перетворюється в нуль. Її перший аргумент – індекс функції Бесселя, другий – індекс кореня. Визначимо функцію mu()
> mu:=n->BesselJZeros(0,n);
BesselJZeros
Власні функції будуть мати наступну структуру.
> R:=(rho,n)->BesselJ(0,rho*mu(n)/L);
BesselJ
Очевидно, що
початковий профіль функції
співпадає з профілем
(оскільки f(0)=0).
> solve(u(rho,0)=v(rho,0)+F(rho)*f(0),v(rho,0));
Оскільки, це значення нульове за умовою задачі, то такою ж має бути умова і для функції T(t), яку визначаємо на даному етапі розв’язування.
> dsolve({diff(T(t),`$`(t,2))=T(t)*a^2*(-lambda^2),T(0)=0}, T(t));
Параметр lambda повинен співпадати з одним із власних значень, які були знайдені вище.
> T:=(t,n)->sin(a*t*mu(n)/L);
Розв’язок для функції v(rho,t) шукаємо у вигляду ряду за власними функціями.
> v:=(rho,t)->Sum(B(n)*R(rho,n)*T(t,n),n=1..infinity);
> v(rho,t);
Коефіцієнти B(n) визначимо, використовуючи початкові умови для похідної по часу від функції v(rho,t).
> F(rho)*D(f)(0);
> Vt0:=-%;
Визначимо процедуру обчислення похідної по часу від функції v(rho,t), яка подана у вигляді ряду за функціями Бесселя. Процедуру визначаємо так, щоб можна було обчислити значення цієї похідної в точці.
> Vt:=proc(rho,t)
> local s;
> diff(v(rho,s),s);
> simplify(subs(s=t,%));
> end proc:
В рамках процедури тимчасова змінна замінюється локальною, за цією змінною обчислюється похідна, потім в одержаному виразі локальна змінна замінюється на аргумент, який вказаний другим під час виклику процедури, і після цього вираз спрощується.
> Vt(rho,t);
> Vt(rho,0);
> eqn:=B(n)*a*BesselJZeros(0,n)=int(Vt0*rho*BesselJ(0,
rho*BesselJZeros(0,n)/L),rho=0..L)/int(rho*BesselJ(0, rho*BesselJZeros(0,n)/L)^2,rho=0..L);
> B:=unapply(solve(eqn,B(n)),n);
> u:=(rho,t)->v(rho,t)+F(rho)*f(t);
Таким чином, розв’язок задачі має наступний вигляд.
> u(rho,t);
Далі присвоюємо числові значення параметрам і створюємо анімаційну картинку.
> A:=1;
> a:=1;
> omega:=2*Pi;
> L:=1;
Однак, необхідно спочатку перевизначити деякі залежності. Так, при обчисленні коефіцієнтів B(n) інтеграл будемо обчислювати наближено, для чого введемо нові коефіцієнти.
> C:=n->evalf(B(n));
Сам нескінченний ряд обриваємо на доданку з номером N і замінюємо коефіцієнти B(n) на C(n).
> U:=(rho,t,N)->sum(C(n)*R(rho,n)*T(t,n),n=1..N)+F(rho)*f(t);
Наприклад, так виглядає наша числова функція при врахуванні тільки трьох доданків в ряді.
> U(rho,t,3);
Крім того, перш ніж будувати зображення, необхідно врахувати що стандартна циліндрична система координат в Maple визначена таким чином, що побудова поверхонь в ній здійснюється, коли задається залежність радіуса від висоти, а не висоти від радіуса, як зазвичай прийнято. Тому визначимо власну систему координат cylind так, щоб залежність задавалась звичним чином.
> addcoords(cylind,[z,r,phi],[r*cos(phi),r*sin(phi),z]);
Тепер створюємо анімацію. Необхідно відзначити, що хоча функція і не залежить від кута phi, діапазон його зміни все ж необхідно задавати. Крім того, важливо на якому проміжку часу відображається функція і скільки при цьому складових доданків ряду. На рисунках можна прослідкувати як змінюється стан системи на протязі 1 секунди (рис. 3.10-3.13).
>plots[animate3d](U(rho,t,5),rho=0..1,phi=0..2*Pi,t=0..1, cords=cylind,axes=FRAME,titlefont=[HELVETICA,BOLD,12]);
Коливання круглої мембрани. Перший кадр (початковий).
Коливання круглої мембрани. Третій кадр.
Коливання круглої мембрани. П’ятий кадр.
Коливання круглої мембрани. Восьмий кадр.
З процедурами відтворення анімації необхідно працювати дуже обережно. Наприклад. Якщо збільшити часовий проміжок, то одержимо дещо іншу ситуацію. На рисунках (3.14–3.16) можна бачити три кадри підряд, на яких зображено стан мембрани в різні моменти часу.
>plots[animate3d](U(rho,t,5),rho=0..1,phi=0..2*Pi,t=0..5, coords=cylind,axes=FRAME,titlefont=[HELVETICA,BOLD,12]);
Коливання круглої мембрани. Другий кадр.
Коливання круглої мембрани. Третій кадр.
Коливання круглої мембрани. Четвертий кадр.