Тема 2 ЗВЕДЕННЯ РІВНЯНЬ ДРУГОГО ПОРЯДКУ ДО КАНОНІЧНОГО ВИГЛЯДУ ЗА ДОПОМОГОЮ ЗАМІНИ ЗМІННИХ 5 страница
Задовольняємо початкові умови
> U1:=U(x,0);
> f:=x->(x<=0,0,x<=Pi/2,x,x<=Pi,Pi-x,o);
Умови Діріхле для функції 
виконуються (функція неперервна і кусково-монотонна на проміжку 
), тому її можна розкласти в ряд Фур’є.
> U(x,o)=f(x);
> assume(n::posint);
> f1(x):=x;
> f2(x):=Pi-x;
> A[n]:=2/Pi*{int(f1(x)*sin(n*x),x=0..Pi/2)+ +int(f2(x)*sin(n*x),x=Pi/2..Pi)};
> simplify(%);
> d:=combine(%,tgig);
> sin(1/2*Pi*n):=(-1)^n;
> C[n]:=d;
> U1;
- розв’язок рівняння коливання струни, який задовольняє початкові і крайові умови
Зауваження: для парних значень 
коефіцієнт 
дорівнює нулю, тому для 
коефіцієнт розкладу набуде вигляду
, відповідно
- розв’язок рівнянняколивання струни, який задовольняє початкові і крайові умови.
Проаналізуємо одержаний розв’язок, відобразивши його графічно. Оскільки, ряд для функції U(x,t) нескінченний, то для графічного відображення необхідно залишити скінченну кількість доданків. Відповідний вираз визначимо наступним чином (коефіцієнти записані в явному вигляді, N – число доданків в ряді).
> S:=(x,t,N)->sum(2*(2*(-1)^n/n^2)*sin(n*x)/Pi,n=1..N);
Тепер за допомогою процедури animate () відтворимо процес коливання струни (рис.3.1).
> plots[animate](S(x,t,10),x=0..15,t=0..1, numpoints=100, titlefont=[HELVETICA,BOLD,12]);
Приклад 3.3 Знайти закон розподілу температури 
всередині стержня, розміщеного на відрізку 
, якщо в початковий момент температура всередині стержня була рівною 
. На кінцях стержня підтримується нульова температура.
Знайти розв’язок рівняння 
, яке задовольняє крайові умови 
, 
і початкову умову 
.
Розв’язування
Рівняння, яке розглядається, є рівнянням гіперболічного типу. Подібні задачі виникають при описування процесу поширення хвиль. Функція u(x,t) визначає відхилення точки струни з координатою х в момент часу t. Нульові початкові умови відповідають ситуації, коли кінці струни закріплені. Відповідне хвильове рівняння (а – параметр задачі) має наступний вигляд.
> Eq:=diff(u(x,t),t)=a^2*diff(u(x,t),x$2);
> Eq;
> a:=1;Eq;
Розв’язок рівняння будемо шукати у вигляді ряду за власними функціями.
> U:=(x,t)->Sum(A[n]*(exp((-1)*a^2*(Pi*n/l)^2*t))*sin(Pi*n*x/l) ,n=1..infinity);
> a:=1;
> U(x,t);
> U1:=U(x,0);
Невідомі коефіцієнти розкладу А[n] знаходяться із початкової умови. Відповідно до умови задачі, в початковий момент профіль струни визначається наступною функцією.
> f:=x->(x<=0,0,x<=l/2,x,x<=l,l-x,o);
> U(x,o)=f(x);
> assume(n::posint);
> f1(x):=x;
> f2(x):=l-x;
> A1[n]:=2/l*{int(f1(x)*sin(Pi*n*x/l),x=0..l/2)};
> A2[n]:=2/l*{int(f2(x)*sin(Pi*n*x/l),x=l/2..l)};
> sin(1/2*Pi*n):=(-1)^n;
> A1[n];A2[n];
> A[n]:=4*l/(Pi^2*n^2)*sin(Pi*n/2);
> U(x,t);
-розв’язок рівняння теплопровідності, який задовольняє початкові і крайові умови.
> S:=(x,l,N)->sum(4*sin(Pi*n*x/l)*l*(-1)^n/Pi^2*n^2,n=1..N);
Оскільки, ряд для функції u(x,t) нескінченний, то його необхідно обмежити – необхідно залишити скінченну кількість доданків.
Тепер за допомогою процедури animate () відтворимо процес коливання струни (рис.3.2-3.6).
> plots[animate](S(x,l,5),x=0..100,l=1..5,view=-100..100, numpoints=100,titlefont=[HELVETICA,BOLD,12]);
Момент часу
t=0 c.
Момент часу
t=1 c.
Момент часу
t=5 c.
Момент часу
t=10 c.
Момент часу
t=14 c.
Приклад 3.4 Розв’язати задачу про коливання прямокутної мембрани:
Розв’язування
>Еqn:=diff(u(x,y,t),t$2)=a^2*(diff(u(x,y,t),x$2)+ +diff(u(x,y,t), y$2));
> f:=(x,y)->A*x*y*(l-x)*(L-y);
> pdsolve(Eqn,HINT=X(x)*Y(y)*T(t));
> dsolve({diff(X(x),`$`(x,2))=-lambda^2*X(x),X(0)=0},X(x));
> _EnvAllSolutions:=true;
> solve(sin(lambda*l)=0,lambda);
> about(_Z1);
Originally _Z1, renamed _Z1~:
is assumed to be: integer
> nu:=n->Pi*n/l;
> mu:=n->Pi*n/L;
> X:=(x,n)->sin(x*nu(n));
> Y:=(y,m)->sin(y*mu(m));
>dsolve({diff(T(t),`$`(t,2))=a^2*_c[l]*T(t)+ +a^2*_c[2]*T(t),D(T)(0)=0},T(t));
> T:=(t,n,m)->cos(a*t*sqrt(nu(n)^2+mu(m)^2));
>S:=(x,y,t,N,M)->Sum(Sum(U(n,m)*X(x,n)*Y(y,m)*T(t,n,m), n=1..N),m=1..M);
> u:=(x,y,t)->S(x,y,t,infinity,infinity);
> u(x,y,t);
> u(x,y,0);
`
> assume(n::posint,m::posint);
>U:=(n,m)->int(int(f(x,y)*X(x,n)*Y(y,m),y=0..L), x=0..l)/int(int(X(x,n)^2*Y(y,m)^2,y=0..L),x=0..l);
> u(x,y,t);
Щоб уявити, як буде коливатися мембрана, створимо анімаційну картинку. Але перш за все присвоїмо числові значення параметрам.
> a:=1;
> l:=1;
> L:=1;
> A:=1;
Створюємо анімацію ( в сумі для функції залишаємо 10 доданків по кожному із індексів сумування – всього 100) за допомогою процедури animate3d трьохмірної графіки із пакета plots (рис.3.7-3.9).
>plots[animate3d](S(x,y,t,10,10),x=0..1,y=0..1,t=0..sqrt(2), axes=FRAME,style=HIDDEN,color=BLACK,orientation=[50,60]);
Коливання квадратної мембрани. Момент часу t=1 c.
Коливання квадратної мембрани. Момент часу t=3 c.
Коливання квадратної мембрани. Момент часу t=4 c.
Приклад 3.5 Розв’язати задачу про коливання круглої мембрани: 
, 
, 
.
Розв’язування
Ця задача, в порівнянні з попередньою, має дві принципові особливості. По перше, шукана функція залежить тільки від двох змінних – часу 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);
