|
Работа опубликована в интернет-журнале "Математика в Вузе", №3
Вернуться на страницу <Методические разработки>
| Задано дифференциальное уравнение с неизвестной функцией y(x) |
|
|
|
(1) |
| где |
|
|
|
|
| и краевые условия (слева - 1 рода, справа - 3 рода) |
|
|
|
(2) |
|
|
(3) |
| где a=0, b=1, r=10. Требуется разработать и опробовать алгоритм приближенного решения таких задач методом конечных элементов. |
|
|
|
|
Построение алгоритма
|
|
|
|
Для приближенного решения уравнения (6) выберем систему базисных функций, линейной комбинацией которых будет искомое приближенное решение. Для этого разобьем промежуток [a,b] на части (вообще, неравные) (узловыми) точками a= x0 , x1 , x2 ,.. xn = b . Кусочно-линейную базисную функцию  при k=1,..(n-1) зададим формулой |
|
|
|
(4) |
| Она имеет нулевое значение везде, кроме промежутка (xk-1 , xk-1), в точке xk она равна 1, график ее имеет вид треугольника (см. ниже). |
(5) |
| Для k=0 и для k=n базисные функции задаются второй или первой формулами из (4). |
|
Таким образом, при k=1,3,..n базисные функции не равны нулю лишь на двух интервалах, соприкасающихся в точке xk , а базисные функции и не равны нулю лишь на первом и последнем интервалах. |
(6) |
| Вместо искомой функции y(x) будем искать лишь ее значения в узловых точках и в связи с этим будем использовать в уравнениях вместо y(x) кусочно-линейную функцию |
|
|
|
(7) |
| где vi - ее значение в узловой точке xi. |
|
| Преобразуем дифференциальное уравнение (1), умножив обе его части на базисную функцию и проинтегрировав обе части равенства в границах промежутка [a,b]: |
|
|
|
(8) |
| Эта операция с точностью до числового множителя совпадает с вычислением ортогональной проекции обеих частей равенства на базисный вектор в пространстве функций, интегрируемых с квадратом на промежутке [a,b] и со скалярным произведением |
|
|
|
(9) |
| Но доставляемые такой операцией числа не совпадают с координатами проектируемой функции в нашем базисе, потому что среди базисных функций ортогональными являются лишь те, номера которых отличаются более, чем на единицу. |
|
| Интегрируя в (08) по частям первое слагаемое, с учетом краевого условия (3) получим |
|
|
|
(10) |
| Отметим, что производная кусочно-линейной базисной функцииявляется кусочно-постоянной с разрывами конечной величины в узловых точках, что не создает проблем для ее интегрирования. |
|
| Подставим сюда вместо y(x) функцию (7) и используем свойство (5) для уменьшения промежутка интегрирования. В силу (7) заменим y(b) на vn . Тогда для k = 0,1,..n получим систему из (n+1) уравнений относительно n+2 неизвестных { y' (a), v0 , v1 , v2 , v3 ,.. vn }. |
|
|
|
(11a) |
|
|
(11b) |
|
|
(11c) |
| Добавив сюда краевое условие (2), котороe при замене функции y(x) на v</<U>7) приобретет вид |
|
|
|
(12) |
| мы получим замкнутую систему n+2 линейных алгебраических уравнений относительно n+2 неизвестных y' (a), v0 , v1 , v2 , v3 ,.. vn |
|
|
|
| Выясним разрешимость полученной системы. Уравнение (12) сразу уменьшает число неизвестных. Уравнение (11a) пригодно для отыскания неизвестной y' (a), поскольку в других уравнениях ее нет. Осталась замкнутая подсистема n уравнений (11b) - (11c) относительно n неизвестных { v1 , v2 , v3 ,.. vn }, имеющая трехдиагональную матрицу коэффициентов A. |
|
|
|
| Обозначим < f | g > интеграл вида |
|
|
|
(13) |
| Очевидно, что это - билинейный симметричный функционал. Благодаря положительности функций p(x) > 0 и q(x) > 0 он еще и невырожденный: |
|
|
|
(14) |
| и поэтому может служить скалярным произведением функций. Тогда коэффициенты уравнений (11b) - (11c) |
|
|
|
(15) |
| По формулам (15) видно, что если бы не было слагаемого rp(b) в коэффициенте Ann , то матрица A была бы матрицей Грама системы линейно независимых векторов . Скалярный квадрат < v | v > кусочно-линейной функции v вида (7) с одной стороны является квадратичной формой относительно величин v0 , v1 , v2 , v3 ,.. vn и матрица коэффициентов этой формы - как раз матрица A без добавки rp(b) в ее элементе Ann. C другой стороны он обладает свойством строгой положительной определенности (14) в пространстве всех функций вида (7) (строгая положительность следует из того факта, что функции вида (7) не могут отличаться от нуля лишь на множестве нулевой меры; скалярный квадрат для функций (7) равен 0 только если функция везде равна 0, т.е. равны 0 все v0 , v1 , v2 , v3 ,.. vn ). |
|
| Получается, что A без добавки rp(b) в ее элементе Ann была бы положительно определенной симметричной матрицей. Добавка rp(b) в элементе Ann означает лишь добавление квадратичного слагаемого rp(b) vn2 к строго положительной квадратичной форме < v | v > . Очевидно, эта добавка из-за положительности коэффициента rp(b) не нарушает строгой положительности квадратичной формы. Не нарушает она и симметрии матрицы. |
|
| В итоге мы делаем вывод: A - симметричная и строго положительно определенная матрица. Из теории квадратичных форм известно, что любая квадратичная форма приводится ортогональным преобразованием коэффициентов к сумме квадратов с коэффициентами, равными собственным числам матрицы. Отсюда мы делаем вывод, что все собственные числа матрицы A положительны, т.е. не равны нулю. Значит, теоретически система уравнений (11b) - (11c) разрешима. |
|
| Ее трехдиагональность позволяет применять метод прогонки. |
|
|
|
|
Тестовый пример
|
|
|
Использован Maple 6 (tm), его команды - красным цветом (так в нем принято).
Программа позволяет заменить базисные функции другими.
При необходимости видеть результат команды ставьте после нее вместо двоеточия точку с запятой
|
|
> restart;with(linalg):with(plots):
Warning, the protected names norm and trace have been redefined and unprotected
Warning, the name changecoords has been redefined
|
|
|
Границы, N (число подинтервалов), коэффициенты уравнений (1, 3)
> a:=0; b:=1; N:=10; p:=1/(1+x^2); q:=1+exp(-x^2); r:=10;
|
|
|
Сами задаем решение уравнения (1)
> y:=(x-a)*exp(c*x);
Подбираем параметр с для выполнимости второго краевого условия (3)
> eq:=subs(x=b,diff(y,x)+r*y): c:=solve(eq=0,c);
> pic1:=plot(y,x=a..b):
|
|
|
Контроль выполнения краевых условий (2, 3)
> subs(x=a,y); subs(x=b,diff(y,x)+r*y);
|
|
|
Вычисляем правую часть уравнения (1)
> f:=diff(-p*diff(y,x),x)+q*y;
|
|
|
Дробление промежутка [a,b] (может быть и неравномерным)
> xk:=array(-1..N+1):
> for k from -1 to N+1 do xk[k]:=a+k*(b-a)/N: end do:
|
|
|
Базисные функции (4) и их графики
> fik:=(k,x)->piecewise(xk[k-1]<x and x<=xk[k],(x-xk[k-1])/(xk[k]-xk[k-1]), xk[k]<x and x<xk[k+1],(x-xk[k+1])/(xk[k]-xk[k+1]),0);
> plot({fik(0,x),fik(1,x),fik(2,x),fik(N,x)},x=a..b);
Рис. 1. Базисные функции (первые три и последняя).
|
|
| Составляем матрицу коэффициентов системы уравнений (11b, 11c) и вектор их правых частей |
|
| > A:=matrix(N,N,0): ff:=vector(N): |
|
| Значения правых частей системы уравнений |
|
|
> for k from 1 to (N-1) do ff[k]:=evalf(Int(f*fik(k,x),x=xk[k-1]..xk[k+1])) end do:
> ff[N]:=evalf(Int(f*fik(N,x),x=xk[N-1]..xk[N])):
|
|
| Матрица A ( 15 ) уравнений ( 11b, 11c) |
|
|
> for k from 1 to (N-1) do dfk1:=diff(fik(k,x),x): dfk2:=diff(fik(k+1,x),x): fk1:=fik(k,x): fk2:=fik(k+1,x): A[k,k]:=evalf(Int(p*dfk1*dfk1+q*fk1*fk1,x=xk[k-1]..xk[k]))+ evalf(Int(p*dfk1*dfk1+q*fk1*fk1,x=xk[k]..xk[k+1])): A[k,k+1]:=evalf(Int(p*dfk1*dfk2+q*fk1*fk2,x=xk[k]..xk[k+1])): A[k+1,k]:=A[k,k+1] end do:
> pb:=subs(x=b,p): dfN:=diff(fik(N,x),x):fkN:=fik(N,x): A[N,N]:=r*pb+evalf(Int(p*dfN^2+q*fkN^2,x=xk[N-1]..xk[N])):
|
|
|
> evalm(A): evalm(ff):
Решение системы A*V=ff
> v:=linsolve(A,ff):
> pic2:=listplot([[a,0],seq([xk[k],v[k]],k=1..N)]):
> display([pic1,pic2]);
|
|
|
|
|
|
Рис 2. Точное (красным цветом) и приближенное решения (N=10)
При более мелком дроблении совпадение более хорошее. Так, например при N=50 получилось:
|
|
|
|
|
|
Рис.3. Точное (красным цветом) и приближенное решения (N=50).
Графики почти сливаются.
Программа позволяет задавать неравномерное дробление промежутка [a,b]. За счет этого можно, получив на крупном дроблении примерный вид решения, измельчить дробление на участках, где решение имеет большую производную. Это позволит обойтись меньшим числом неизвестных vi.
|
|
Наверх
Вернуться на страницу <Методические разработки>
|
