Тема 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);