Решение уравнения пуассона в mathcad

Глава 5. Решение дифференциальных уравнений

5.9 Решение эллиптических уравнений (Лапласа и Пуассона)

С помощью двух встроенных функций, multigrid и relax , можно решить простейшие случаи уравнения Пуассона. Это уравнение очень часто используется в технике, например, для описания полей напряжений и деформаций в плоской задаче теории упругости, в задачах теплопроводности, гидроаэродинамики, электростатики. В MathCAD для численного решения уравнения Пуассона используется метод конечных разностей.

Уравнение Пуассона имеет вид

Если правая часть уравнения равна нулю, такое уравнение называется уравнением Лапласа

На квадратной области уравнение Пуассона представляется в виде

Численное решение ищется в MathCAD только на квадратной области, состоящей из ( n +1) ( n +1) точек. Поэтому граничные условия должны быть определены пользователем для всех четырех сторон квадрата. Самый простой ( и наиболее часто используемый) вариант – это нулевые граничные условия (уравнение Лапласа). В таком случае можно использовать функцию multigrid .

Обращение к функции:

Multigrid ( M , n cycle ),

где M – квадратная матрица размером 1 2 n , которая содержит значения правой части уравнения Пуассона в соответствующей точке квадратной области; n cycle – число циклов на каждом уровне итерации функции multigrid . Значение n cycle =2 обычно дает хорошую аппроксимацию решения.

Пример использования функции multigrid приведен на рис. 5.21.

Обнуление предыдущих значений М

Обнуление матрицы правых частей уравнения

Три точечных источника

Значения правой части уравнения Пуассона

Рис. 5. 21 Решения уравнения Лапласа с помощью функции multigrid

Если граничные условия по сторонам квадрата ненулевые, необходимо использовать функцию relax .

Обращение к функции:

relax(a, b, c, d, f, u, rjac),

где a , b , c , d , e – квадратные матрицы одинакового размера, содержащие коэффициенты аппроксимирующего уравнения; f – квадратная матрица, содержащая значения правой части уравнения в каждой точке области, где ищется решение; u – квадратная матрица, содержащая граничные значения решения на границе квадратной области и начальное приближение для решения внутри области; rjac – спектральный радиус итераций Якоби. Это число между 0 и 1, которое управляет сходимостью процесса релаксации.

Использование этой функции требует глубокого знания метода конечных разностей для составления указанных матриц. Пример использования функции relax приведен на рис. 5.22 и 5.23.

Определим размеры сетки

Введем 5 квадратных матриц для коэффициентов a, b, c, d, e ,

входящих в сеточную аппроксимацию уравнения Пуассона

чем больше эти коэффициенты, тем меньше шаг решения,

тем точнее результат. Можно взять

Задает положение и интенсивность источника

Пусть во всех узлах значения правой части одинаковы

Граничные условия

на верхней границе

на нижней границе

по бокам

Поменяйте условия. Включите серые выражения.

Это аналог функции multigrid

Спектральный радиус Якоби r

Решение уравнения Пуассона

Рис. 5. 22 Решение уравнения Пуассона с помощью функции relax

Рис. 5. 23 Результаты решение уравнения Пуассона с помощью функции relax

Решение уравнения пуассона в mathcad

11.3.2. Эллиптические уравнения

Решение эллиптических уравнений в частных производных реализовано только для единственного типа задач — двумерного уравнения Пуассона.

Это уравнение содержит вторые производные функции u(х,у) по двум пространственным переменным;

Уравнение Пуассона описывает, например, распределение электростатического поля u(х,у) в двумерной области с плотностью заряда f (x,y) , или (см. разд. 11.1.2) стационарное распределение температуры u(х,у ) на плоскости, в которой имеются источники (или поглотители) тепла с интенсивностью f (х,у) .

Несмотря на то, что применение встроенных функций, описанных в данном разделе, анонсировано разработчиками Mathcad только для уравнения Пуассона, их можно применять и для решения других уравнений, даже необязательно эллиптического типа. О том, как осуществить такие расчеты, написано в конце данного раздела.

Уравнение Пуассона с нулевыми граничными условиями

Корректная постановка краевой задачи для уравнения Пуассона требует задания граничных условий. В Mathcad решение ищется на плоской квадратной области, состоящей из (M+1)х(M+1) точек. Поэтому граничные условия должны быть определены пользователем для всех четырех сторон упомянутого квадрата. Самый простой вариант — это нулевые граничные условия, т. е. постоянная температура по всему периметру расчетной области. В таком случае можно использовать встроенную функцию multigrid :

multigrid (F,ncycie) — матрица решения уравнения Пуассона размера (м+1)х(м+1) на квадратной области с нулевыми граничными условиями:

F — матрица размера (M+1)х(M+1) , задающая правую часть уравнения Пуассона;
ncycle — параметр численного алгоритма (количество циклов в пределах каждой итерации).

Сторона квадрата расчетной области должна включать точно N=2n шагов, т. е. 2n+1 узлов, где n — целое число.

Параметр численного метода ncycle в большинстве случаев достаточно взять равным 2. Листинг 11.7 содержит пример использования функции multigrid для расчета краевой задачи на области 33х33 точки и точечным источником тепла в месте, задаваемом координатами (15,20) внутри этой области.

Листинг 11.7. Решение уравнения Пуассона с нулевыми граничными условиями

В первой строке листинга задается значение м=32, в двух следующих строках создается матрица правой части уравнения Пуассона, состоящая из всех нулевых элементов, за исключением одного, задающего расположение источника. В последней строке матрице с присваивается результат действия функции multigrid . Обратите внимание, первый ее аргумент сопровождается знаком «минус», что соответствует записи правой части уравнения Пуассона (11.11). Графики решения показаны на рис. 11.17 и 11.18 в виде трехмерной поверхности и линий уровня соответственно.

Рис. 11.17. График поверхности решения уравнения Пуассона (продолжение листинга 11.7)

Рис. 11.18. График линий уровня решения уравнения Пуассона (продолжение листинга 11.7)

Уравнение Пуассона с произвольными граничными условиями

В более сложных случаях, например, для решения краевой задачи с ненулевыми условиями на границах, следует использовать другую встроенную функцию relax , имеющуюся в Mathcad:

relax(a,b,c,d,e,F,v,rjac) — матрица решения дифференциального уравнения в частных производных на квадратной области, полученного с помощью алгоритма релаксации для метода сеток:

a,b,c,d,e — квадратные матрицы коэффициентов разностной схемы, аппроксимирующей дифференциальное уравнение;
F — квадратная матрица, задающая правую часть дифференциального уравнения;
v — квадратная матрица граничных условий и начального приближения к решению;
rjac — параметр численного алгоритма (спектральный радиус итераций Якоби);

Суть алгоритма релаксации сводится к тому, что в ходе итераций происходит проверка уравнений и соответствующая коррекция значений искомой функции в каждой точке. Если начальное приближение выбрано удачно, то можно надеяться, что алгоритм сойдется («релаксирует») к правильному решению.

Все матрицы, задающие как коэффициенты разностной схемы а, b, с, d, e , граничные условия v, так и само решение F , должны иметь одинаковый размер (M+1)х(M+1) , соответствующий размеру расчетной области. При этом целое число м обязательно должно быть степенью двойки: м=2п.

Решение уравнения Пуассона с тремя источниками разной интенсивности при помощи функции relax приведено в листинге 11.8.
Листинг 11.8. Решение уравнения Пуассона с помощью функции relax

Первые три строки имеют тот же смысл, что и в предыдущем листинге. Только вместо одного источника тепла взято их другое распределение — один сильный источник, один более слабый и один сток тепла. В следующих шести строках задаются коэффициенты разностной схемы. Отложим их обсуждение до последнего раздела этой главы, ограничившись утверждением, что для решения уравнения Пуассона коэффициенты должны быть взяты именно такими, как показано в листинге 11.8. В предпоследней строке задана матрица нулевых граничных условий и нулевых начальных приближений, а в последней матрице с присваивается результат действия функции relax . График полученного решения в виде линий уровня показан на рис. 11.19.

Рис. 11.19. Решение уравнения Пуассона с помощью функции relax (продолжение листинга 11.8)

Разностная схема для решения уравнения Пуассона

Несмотря на отсутствие сведений в справочной системе Mathcad о решении других линейных дифференциальных уравнений в частных производных, кроме уравнения Пуассона, сделать это возможно с помощью той же функции relax (см. предыдущий разд.). Для этого нужно правильным образом задать коэффициенты разностной схемы.

Начнем с пояснения выбора этих коэффициентов (см. листинг 11.8) для уравнения Пуассона. Согласно основным идеям метода сеток (см. разд. 11.2), для дискретизации обеих пространственных производных в уравнении (11.12) следует использовать по три соседних узла вдоль каждой из координат. Поэтому уравнение Пуассона (11.12) может быть записано в разностной форме при помощи шаблона типа «крест» (рис. 11.20). В этом случае, после приведения подобных слагаемых в разностных уравнениях коэффициенты разностной схемы будут такими, как показано возле узлов шаблона на этом рисунке (аналогичные коэффициенты для явной и неявных схем решения уравнения теплопроводности см. на рис. 11.6 и 11.11 соответственно).

Теперь, если вы сравните полученные числа с константами, которые присвоены элементам матриц-аргументов функции relax (см. листинг 11.8), то увидите, что они как раз и описывают вычисленные нами только что коэффициенты разностной схемы «крест». Таким образом, нетрудно сообразить, что с помощью встроенной функции relax можно решать и другие линейные дифференциальные уравнения в частных производных, которые можно аппроксимировать схемой типа «крест» или схемой, являющейся ее составной частью. Конечно, для того чтобы использовать эту встроенную функцию для другого уравнения, необходимо будет составить соответствующую разностную схему.

Рис. 11.20. Шаблон аппроксимации уравнения Пуассона «крест»

Решение уравнения диффузии тепла при помощи функции relax

Приведем пример применения встроенной функции relax для решения другого уравнения в частных производных (т. е. не уравнения Пуассона, для которого она изначально предназначена). Вычислим при помощи этой функции решение уже хорошо нам знакомого однородного линейного уравнения теплопроводности (см. разд. 11.1.2). Будем использовать явную разностную схему, шаблон которой изображен на рис. 11.6. Для того чтобы «приспособить» для явной схемы функцию relax , требуется только задать ее аргументы в соответствии с коэффициентами, показанными на шаблоне (см. рис. 11.6). Программа, реализующая таким способом явную схему, представлена на листинге 11.9. Число Куранта в этом листинге обозначено переменной с, как и положено явной разностной схеме, она выдает устойчивое решение только для C .

Листинг 11.9. Решение уравнения теплопроводности при помощи функции relax

Результат действия программы листинга 11.9 показан на рис. 11.21 в виде трехмерной поверхности. Если сравнить рис. 11.21 с рис. 11.4, полученным при расчетах по запрограммированной разностной схеме, то в графиках рис. 11.4 нетрудно узнать сечения этой поверхности плоскостями t=const . Еще раз подчеркнем, что использовать встроенную функцию можно только для тех уравнений, которые допускают построение разностной схемы типа «крест» (см. рис. 11.17) или составного фрагмента этой схемы.

Рис. 11.21. Решение уравнения теплопроводности с помощью функции relax (продолжение листинга 11.9)

Эллиптические уравнения

Решение эллиптических уравнений в частных производных реализовано только для единственного типа задач — двумерного уравнения Пуассона. Это уравнение содержит вторые производные функции u(х,у) по двум пространственным переменным:

Уравнение Пуассона описывает, например, распределение электростатического поля u(х,у) в двумерной области с плотностью заряда f (х,у) или (см. разд. 13.1.2) стационарное распределение температуры u(х,у) на плоскости, в которой имеются источники (или поглотители) тепла с интенсивностью f (х,у).

Уравнение Пуассона с нулевыми граничными условиями

Корректная постановка краевой задачи для уравнения Пуассона требует задания граничных условий. В Mathcad решение ищется на плоской квадратной области, состоящей из (м+1)х(м+1) точек. Поэтому граничные условия должны быть определены пользователем для всех четырех сторон упомянутого квадрата. Самый простой вариант — это нулевые граничные условия, т. е. постоянная температура по всему периметру расчетной области. В таком случае можно использовать встроенную функцию multigrid.

muitigrid(F,ncycle) — матрица решения уравнения Пуассона размера (M+1)х(M+1) на квадратной области с нулевыми граничными условиями;
- F — матрица размера (M+1)X(M+1), задающая правую часть уравнения Пуассона;
- ncycie — параметр численного алгоритма (количество циклов в пределах каждой итерации).

Сторона квадрата расчетной области должна включать точно M=2 n шагов, т.е. 2 n +1 узлов, где n — целое число.

Параметр численного метода ncycie в большинстве случаев достаточно взять равным 2. Листинг 13.6 содержит пример использования функции multigrid для расчета краевой задачи на области ззхзз точки и точечным источником тепла в месте, задаваемом координатами (15,20) внутри этой области.

Листинг 13.6. Решение уравнения Пуассона с нулевыми граничными условиями.

В первой строке листинга задается значение M=32, в двух следующих строках создается матрица правой части уравнения Пуассона, состоящая из всех нулевых элементов, за исключением одного, задающего расположение источника. В последней строке матрице G присваивается результат действия функции multigrid. Обратите внимание, первый ее аргумент сопровождается знаком «минус», что соответствует записи правой части уравнения Пуассона (11). Графики решения показаны на рис. 13.16 и 13.17 в виде трехмерной поверхности и линий уровня, соответственно.

Уравнение Пуассона с произвольными граничными условиями

В более сложных случаях, например для решения краевой задачи с ненулевыми условиями на границах, следует использовать другую встроенную функцию relax, имеющуюся в Mathcad.

reiax(a,b,c,d,e,F,v,rjac) — матрица решения дифференциального уравнения в частных производных на квадратной области, полученного с помощью алгоритма релаксации для метода сеток;
- a,b,c,d,e— квадратные матрицы коэффициентов разностной схемы, аппроксимирующей дифференциальное уравнение;
- F — квадратная матрица, задающая правую часть дифференциального уравнения;
- v — квадратная матрица граничных условий и начального приближения к решению;
- rjac — параметр численного алгоритма (спектральный радиус итераций Якоби).

Рис. 13.16. График поверхности решения уравнения Пуассона (листинг 13.6)

Рис. 13.17. График линий уровня решения уравнения Пуассона (листинг 13.6)

Параметр численного алгоритма характеризует скорость сходимости итераций. Он должен быть числом от о до 1. В матрице граничных условий v необходимо задать только граничные элементы, исходя из значения краевых условий по периметру расчетной области. Прочие (внутренние) элементы этой матрицы служат для задания начального приближения к решению. Суть алгоритма релаксации сводится к тому, что в ходе итераций происходит проверка уравнений и соответствующая коррекция значений искомой функции в каждой точке. Если начальное приближение выбрано удачно, то можно надеяться, что алгоритм сойдется («срелаксирует») к правильному решению.

Все матрицы, задающие как коэффициенты разностной схемы a,b,c,d,e, граничные условия v, так и само решение F, должны иметь одинаковый размер (M+1)х(M+1), соответствующий размеру расчетной области. При этом целое число м обязательно должно быть степенью двойки: M=2″.

Решение уравнения Пуассона с тремя источниками разной интенсивности при помощью функции relax приведено в листинге 13.7.

Листинг 13.7. Решение уравнения Пуассона с помощью функции relax

Первые три строки имеют тот же смысл, что и в предыдущем листинге. Только вместо одного источника тепла взято их другое распределение — один сильный источник, один более слабый и один сток тепла. В следующих шести строках задаются коэффициенты разностной схемы. Отложим их обсуждение до последнего раздела этой главы, ограничившись утверждением, что для решения уравнения Пуассона коэффициенты должны быть взяты именно такими, как показано в листинге 13.7. В предпоследней строке задана матрица нулевых граничных условий и нулевых начальных приближений, а в последней матрице G присваивается результат действия функции relax. График полученного решения в виде линий уровня показан на рис. 13.18.

Рис. 13.18. Решение уравнения Пуассона с помощью функции relax (листинг 13.7)

Разностная схема для решения уравнения Пуассона

Начнем с пояснения выбора этих коэффициентов (см. листинг 13.7) для уравнения Пуассона. Согласно основным идеям метода сеток (см. разд. «Разностные схемы» этой главы), для дискретизации обеих пространственных производных в уравнении (12) следует использовать по три соседних узла вдоль каждой из координат. Поэтому уравнение Пуассона (12) может быть записано в разностной форме при помощи шаблона типа «крест» (рис. 13.19). В этом случае после приведения подобных слагаемых в разностных уравнениях коэффициенты разностной схемы будут такими, как показано возле узлов шаблона на этом рисунке (аналогичные коэффициенты для явной и неявных схем решения уравнения теплопроводности см. на рис. 13.6 и 13.11, соответственно).

Теперь если Вы сравните полученные числа с константами, которые присвоены элементам матриц-аргументов функции relax (см. листинг 13.7), то увидите, что они как раз и описывают вычисленные нами только что коэффициенты разностной схемы «крест». Таким образом, нетрудно сообразить, что с помощью встроенной функции relax можно решать и другие линейные дифференциальные уравнения в частных производных, которые можно аппроксимировать схемой типа «крест» или схемой, являющейся ее составной частью. Конечно, для того чтобы использовать эту встроенную функцию для другого уравнения, необходимо будет составить соответствующую разностную схему.

Рис. 13.19. Шаблон аппроксимации уравнения Пуассона «крест»

Решение уравнения диффузии тепла при помощи функции relax

источники:

http://sistemair.ru/dok/mathcad12/Glava_11/Index11.htm

http://www.bourabai.ru/einf/mathcad/ch13/index9.html