Методическое пособие 521
.pdfВ результате функционал превратился в функцию трёх переменных 1, 2 и . Необходимое условие минимума этой функции
F 1 0, |
F 2 0 |
с учётом изопараметрического условия (4.10) приводит к системе алгебраических уравнений
> eq:={diff(P,a1)=0,diff(P,a2)=0,int(z(x)^2, x=0..1)=1};
Эта система нелинейная. Попытка решить её аналитическими методами приводит к успеху:
> solve(eq,{a1,a2,lambda});
100
Получилось четыре решения. На самом деле здесь различных решений только два, поскольку наборы 1, 2, которые соответствуют одному и тому же собственному числу , отличаются лишь знаками и тем самым дают противоположные по знаку собственные функции. Такие функции считаются неразличимыми, так как собственные функции задачи Штурма-Лиу- вилля определяются с точностью до постоянного множителя.
Покажем, что тот же результат можно получить, не прибегая к решению нелинейной системы. Для этого при составлении системы не будем включать в неё интегральное условие (4.10).
> eq:=[diff(P,a1)=0,diff(P,a2)=0];
Заметим, что полученная система имеет вид (A–B)u = 0, что соответствует известной в линейной алгебре обобщенной задаче на собственные числа и собственные векторы матрицы A. Такая задача сводится к обычной задаче Wu= u определения собственных чисел и векторов для матрицы W=B–1A, методы решения которой для любого порядка матрицы достаточно хорошо разработаны. Например, широко применяются методы вращений (Якоби), итераций, Ланцоша, степенной метод, QR-алгоритм и др.
Чтобы воспользоваться стандартными возможностями системы Maple для нахождения собственных чисел и векторов, необходимо по заданной системе уравнений явно выписать матрицы A и B. Для этого воспользуемся функцией genmatrix из пакета расширения linalg.
>with(linalg):
>A:=genmatrix(subs(lambda=0,eq),[a1,a2]);
101
> B:=evalm(-genmatrix(diff(eq,lambda),[a1,a2]));
Обратите внимание, как здесь из системы (A – B)u = 0 выделяется нужная матрица: в одном случае предварительно подставляется =0, а в другом – берется частная производная по .
Теперь вычисляем собственные числа обобщенной задачи:
> evalf(eigenvals(A,B));
Для нахождения собственных векторов требуется матри-
ца B–1A:
>G:=evalm(inverse(B)&*A);
>v:=evalf(eigenvectors(G));
v:= [41.99999994, 1., {[-.4999999994, 1.]}], [10.00000000, 1., { [1.,
0]}] |
|
~ |
|
~ |
10 |
42 |
|
Итак, получили два собственных числа 1 |
и 2 |
||
. Им соответствуют собственные векторы |
б1 1;0 T и |
б2 0,5;0 T (в разных версиях Maple значения компонент могут отличаться от приведённых некоторым множителем).
Построим собственную функцию, соответствующую первому собственному числу. В приведенном списке v этому числу соответствует второй составной элемент [10.00000000, 1., {[1., 0]}], в котором 10.00000000 – это собственное число, 1.
– его кратность, {[1., 0]} – собственный вектор. Компоненты последнего и есть искомые коэффициенты 1, 2 аппроксимации (*). Извлечь найденные числовые значения 1, 2 из списка можно, например, так:
>num:=2;# номер элемента списка v, соответствующего 1 10
>pr:=a1=op(v[num,3])[1],a2=op(v[num,3])[2];
102
Внимание: так как функция Maple eigenvectors() выводит собственные числа и соответствующие им собственные векторы в случайном порядке, то при практической реализации переменная num может иметь значение, отличное от приведенного, в частности, равное 1. Поэтому необходимо тщательно отслеживать порядок вывода собственных чисел указанной функцией.
Далее делаем подстановку и тут же нормируем, т.е. де-
|
|
|
|
1 2 |
|
1 |
|
|
|
лим на |
|
~2 |
dx |
|
|
y |
|
||
|
0 |
|
|
> z1:=subs(pr,z(x))/sqrt(int(subs(pr,z(x))^2, x=0..1));
z1 := 5.477225576 x (x–1)
Получили приближение первой собственной функции ~y1(x) 5.4772 x (x 1). Интересно сравнить его с точным ре-
шением y1 2sin x ; для этого строим графики и вычисляем
норму погрешности аппроксимации.
> plot({z1,sqrt(2)*sin(Pi*x)},x=0..1);
Рис. 4.2
103
Судя по характерной «картинке», данные функции находятся в противофазе. Как известно, функции с одинаковой нормой могут отличаться знаками, в данном случае так и произошло, поэтому перед z1 нужно взять «минус». Чтобы автоматизировать этот процесс, можно проверять производную в нуле: если она отрицательна, знак собственной функции сле-
|
|
|
0. |
дует поменять, поскольку ( 2sin x) |
|||
|
|
|
x 0 |
|
|
|
> z1:=`if`(subs(x=0,diff(z1,x))>0,z1,-z1); z1 := –5.477225576 x (x–1)
> plot({z1,sqrt(2)*sin(Pi*x)},x=0..1);
Рис. 4.3
Норму погрешности решения определяем по формуле
|
|
|
|
|
|
|
|
|
1 2 |
|
~ |
|
|
|
1 |
|
|
|
|
1 |
y1 |
|
|
~ |
(x) y1 |
2 |
dx |
|
|
y1 |
|
(y1 |
(x)) |
. |
|||||
|
|
|
|
0 |
|
|
|
|
> delta1:=sqrt(int((z1-sqrt(2)*sin(Pi*x))^2,x=0..1));
1:= .03801985008
104
Аналогично находим вторую собственную функцию, со-
~
ответствующую собственному числу 2 42.
> num:=1; # номер элемента списка v, соответствующего
2 42; может не совпадать с указанным
>pr:=a1=op(v[num,3])[1],a2=op(v[num,3])[2];
>z2:=subs(pr,z(x))/sqrt(int(subs(pr,z(x))^2, x=0..1));
>z2:=`if`(subs(x=0,diff(z1,x))>0,z2,-z2);
z2 := 28.98275350 (-.4999999994+1.x) x (x–1)
> plot({z2,sqrt(2)*sin(Pi*x)},x=0..1);
Рис. 4.4
> delta2:=sqrt(int((z2-sqrt(2)*sin(2*Pi*x))^2, x=0..1));
2:= .1308460450
Таким образом, вторая собственная функция аппрокси-
мируется выражением
~y2(x) 28,983x (x 0,5) (x 1).
Для проверки ортогональности вычислим скалярное произведение (~y1,~y2 ); оно должно быть близко к нулю.
105
> int(z1*z2,x=0..1);
-.3174901574 10-8
Проверка сходимости. С увеличением числа m параметров аппроксимации точность решения задачи должна возрастать. Для проведения такого анализа следует решить несколько задач при различных значениях m и показать, что
~2
k m (k ) , k m 0, k = 1, 2, … ,
где k ~yk yk – норма погрешности аппроксимации собст-
венной функции yk . Для этого в приведенной программе следует модифицировать некоторые строки: там, где определяется вид аппроксимации, система уравнений, список неизвестных и т.п. Например, при m=4 имеем (измененные строки выделены)
>restart;
>with(linalg):
>Digits:=20:
>z:=x->(a1+a2*x+a3*x^2+a4*x^3)*x*(x-1):
>F:=y->int((D(y)(x))^2-lambda*y(x)^2,x=0..1):
>P:=evalf(F(z)):
>eq1:=[diff(P,a1),diff(P,a2),diff(P,a3),diff(P,a4)]:
>A:=genmatrix(subs(lambda=0,eq1),[a1,a2,a3,a4]): B:=-genmatrix(diff(eq1,lambda),[a1,a2,a3,a4]): G:=evalm(inverse(B)&*A): evalf(eigenvals(A,B)): v:=evalf(eigenvectors(G));
v := [39.501552784298, 1., { [-.931856758802112e-1, -.212570779438151, 1.19682639311396, -.797884261947117]}], [102.130250362174, 1., { [-.262091514193518, 1.21430684794896, -1.21430684764023, -.283518439568072e-9]}], [200.498447161152, 1., { [.266620374018506, -2.07461390966088, 4.62411948487108, -3.08274632323329]}],
[9.8697496216058, 1., { [.661684886859476, .749782041997382, -.749782042221858, .166378860069766e-9] }]
106
> num:=4: # номер элемента списка v, соответствующего наименьшему собственному числу; при конкретной реализации может быть другим (1, 2 или 3)
pr:=a1=op(v[num,3])[1],a2=op(v[num,3])[2], a3=op(v[num,3])[3], a4=op(v[num,3])[4]:
>z1:=subs(pr,z(x))/sqrt(int(subs(pr,z(x))^2, x=0..1)):
>z1:=`if`(subs(x=0,diff(z1,x))>0,z1,-z1);
>plot({z1,sqrt(2)*sin(Pi*x)},x=0..1);
>delta1:=sqrt(int((z1-sqrt(2)*sin(Pi*x))^2, x=0..1));
z1 := 6.65573235671404(.661684886859476 + .749782041997382 x –
.749782042221858 x2 + .166378860069766e-9 x3) x (x – 1)
>num:=1: # номер элемента списка v, соответствующего второму по величине собственному числу; при конкретной реализации может быть целым числом от 1 до 4
pr:=a1=op(v[num,3])[1],a2=op(v[num,3])[2], a3=op(v[num,3])[3], a4=op(v[num,3])[4]:
>z2:=subs(pr,z(x))/sqrt(int(subs(pr,z(x))^2, x=0..1)):
>z2:=`if`(subs(x=0,diff(z1,x))>0,z2,-z2);
>plot({z2,sqrt(2)*sin(2*Pi*x)},x=0..1);
>delta2:=sqrt(int((z2-sqrt(2)*sin(2*Pi*x))^2, x=0..1));
z2 := 90.0479177822468 (-.931856758802112e-1 – .212570779438151 x
+1.19682639311396 x2 –.797884261947117 x3) x (x – 1)
Видно, что ошибка вычисления первых двух собственных чисел и собственных функций стала меньше, чем при m = 2.
Действуя аналогично, в данном случае с m = 4 можно также получить приближения двух следующих собственных чисел 3 и 4 , точные значения которых равны 9 2 и 16 2 со-
ответственно, и им отвечающие собственные функции. Очевидно, их точность вычисления будет ниже, чем 1 и 2, y1 и y2. Во-
107
обще говоря, с увеличением номера собственного числа (собственной функции) задачи Штурма-Лиувилля скорость сходимости решения метода Ритца уменьшается. Ниже приводится итоговая таблица, с помощью которой удобно вести контроль сходимости.
|
m = 2 |
m = 3 |
m = 4 |
m = 5 |
m = 6 |
точн. знач. |
1 |
10 |
9.86975 |
9.86975 |
9.869604435 |
9.869604435 |
9.869604401 |
|
0.038020 |
0.0008307 |
0.0008307 |
0.8980 10–5 |
0.8980 10–5 |
|
1 |
|
|
|
|
|
|
2 |
42 |
42 |
39.502 |
39.502 |
39.4784734 |
39.478418 |
|
0.13085 |
0.13085 |
0.009748 |
0.009748 |
0.3396 10–3 |
|
2 |
|
|
|
|
|
|
3 |
|
102.130 |
102.130 |
89.174 |
89.174 |
88.8264 |
3 |
|
0.24808 |
0.24808 |
0.036016 |
0.036016 |
|
4 |
|
|
200.5 |
200.5 |
159.99 |
157.91 |
4 |
|
|
0.37680 |
0.37680 |
0.084127 |
|
Из таблицы видна закономерность, при которой результаты расчетов для отдельной пары k и yk меняются в зависимости от m с шагом 2. Например, собственное значение 1 и погрешность 1 остаются неизменными при m = 3 4 и при m = 5 6; 2 и 2 одинаковые для m = 2 и m=3 и для m = 4 и m = 5 и т.д. Это обстоятельство при анализе сходимости позволяет сократить объем вычислений, выполняя их лишь при четных, либо нечетных номерах m.
В задаче Штурма-Лиувилля могут ставиться граничные условия вида
|
(4.11) |
y (a) 0 y(a) 0, |
|
|
|
y (b) 1y(b) 0, |
|
соответствующие однородным условиям 3-го |
рода, а при |
0 = 1 =0 – однородным условиям 2-го рода. В этом случае вариационная формулировка приводит к минимизации функционала с дополнительными граничными слагаемыми
b
F[y(x)] (p(x)y 2 q(x)y2 y2)dx 1p(b)y2(b) 0p(a)y2(a). (4.12)
a
108
Таким образом, задача сводится к решению задачи Больца, в которой ограничения (4.11) выступают как естественные граничные условия – условия трансверсальности. Аппроксимация при этом не должна заранее удовлетворять какому-либо условию при x = a и x = b. С этой целью обычно берут приближение в виде многочлена
~ |
2x 3x |
2 |
... mx |
m 1 |
. |
y(x) 1 |
|
|
Если на одном конце отрезка a x b задано граничное условие 1-го рода, а на другом – 2-3-рода, то в функционале (4.12) следует убрать одно из двух неинтегральных слагаемых, а от аппроксимации ~y(x) потребовать, чтобы выполнялось однородное условие 1-го рода при любых значениях параметров { i}. В частности, для задачи
|
|
|
|
|
|
|
( p(x) y ) q(x)y y; |
|
|
||||
y(a) 0 |
|
|
|
|
, |
|
, y (b) 1y(b) 0 |
||||||
используется функционал |
|
|
|
|
|
|
b |
|
|
|
|
|
|
F[y(x)] (p(x)y 2 q(x)y2 y2)dx 1 p(b)y2(b), |
||||||
a |
|
|
|
|
|
|
аппроксимация функции берется в виде |
|
|
|
|||
~ |
3x |
2 |
... mx |
m 1 |
)(x a), |
|
y(x) ( 1 2x |
|
|
|
|||
а для задачи |
|
|
|
|
|
|
|
|
|
|
|
|
|
( p(x) y ) q(x)y y; |
|
|
|
|
|
|
y(b) 0, |
||
|
y (a) 0 y(a) 0, |
|||||
таковыми соответственно являются |
|
|
|
|||
b |
|
|
|
|
|
|
F[y(x)] (p(x)y 2 q(x)y2 y2)dx 0 p(a)y2(a), |
||||||
a |
|
|
|
|
|
|
~ |
( 1 2x 3x |
2 |
... mx |
m 1 |
)(x b) . |
|
y(x) |
|
|
Пример 3. Найти первые два собственных значения и соответствующие собственные функции задачи
y |
|
y ; |
y(0) = 0, |
|
|
y (1) y(1) 0. |
109