pdenonlin

(Не рекомендуемый), Решают нелинейную эллиптическую задачу УЧП

pdenonlin не рекомендуется. Использование solvepde вместо этого.

Описание

пример

u = pdenonlin(model,c,a,f) решает нелинейный УЧП

(cu)+au=f

с геометрией граничные условия и конечный элемент сцепляются в model, и коэффициенты cA, и f. В этом контексте, “нелинейном”, означает некоторый коэффициент в cA, или f зависит от решения u или его градиент. Если УЧП является системой уравнений (model.PDESystemSize > 1), затем pdenonlin решает систему уравнений

(cu)+au=f

пример

u = pdenonlin(b,p,e,t,c,a,f) решает УЧП с граничными условиями b, и mesh конечного элемента (pET).

пример

u = pdenonlin(___,Name,Value), для любых предыдущих аргументов, изменяет процесс решения с NameЗначение пары.

[u,res] = pdenonlin(___) также возвращает норму остаточных значений шага Ньютона res.

Примеры

свернуть все

Решите минимальную поверхностную задачу. Поскольку эта проблема имеет нелинейный c коэффициент, используйте pdenonlin решить его.

Создайте модель и включайте круговую геометрию с помощью встроенного circleg функция.

model = createpde;
geometryFromEdges(model,@circleg);

Установите коэффициенты.

a = 0;
f = 0;
c = '1./sqrt(1+ux.^2+uy.^2)';

Установите граничное условие Дирихле со значением x2.

boundaryfun = @(region,state)region.x.^2;
applyBoundaryCondition(model,'edge',1:model.Geometry.NumEdges,...
                      'u',boundaryfun,'Vectorized','on');

Сгенерируйте mesh и решите задачу.

generateMesh(model,'GeometricOrder','linear','Hmax',0.1);
u = pdenonlin(model,c,a,f);
pdeplot(model,'XYData',u,'ZData',u)

Figure contains an axes. The axes contains an object of type patch.

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

Создайте геометрию с помощью встроенного circleg функция. Постройте геометрию, чтобы видеть метки ребра.

g = @circleg;
pdegplot(g,'EdgeLabels','on')
axis equal

Figure contains an axes. The axes contains 5 objects of type line, text.

Создайте граничные условия Дирихле со значением x2.Create следующий файл и сохраняют его на вашем пути MATLAB.

function [qmatrix,gmatrix,hmatrix,rmatrix] = pdex2bound(p,e,u,time)
ne = size(e,2); % number of edges
qmatrix = zeros(1,ne);
gmatrix = qmatrix;
hmatrix = zeros(1,2*ne);
rmatrix = hmatrix;
for k = 1:ne
    x1 = p(1,e(1,k)); % x at first point in segment
    x2 = p(1,e(2,k)); % x at second point in segment
    xm = (x1 + x2)/2; % x at segment midpoint
    y1 = p(2,e(1,k)); % y at first point in segment
    y2 = p(2,e(2,k)); % y at second point in segment
    ym = (y1 + y2)/2; % y at segment midpoint
    switch e(5,k)
        case {1,2,3,4}
            hmatrix(k) = 1;
            hmatrix(k+ne) = 1;
            rmatrix(k) = x1^2;
            rmatrix(k+ne) = x2^2;
    end
end

Установите коэффициенты и граничные условия.

a = 0;
f = 0;
c = '1./sqrt(1+ux.^2+uy.^2)';
b = @pdex2bound;

Сгенерируйте mesh и решите задачу.

[p,e,t] = initmesh(g,'Hmax',0.1);
u = pdenonlin(b,p,e,t,c,a,f);
pdeplot(p,e,t,'XYData',u,'ZData',u)

Figure contains an axes. The axes contains an object of type patch.

Решите нелинейную 3-D задачу с нетривиальной геометрией.

Импортируйте геометрию из BracketWithHole.stl файл. Постройте геометрию и столкнитесь с метками.

model = createpde();
importGeometry(model,'BracketWithHole.stl');
figure
pdegplot(model,'FaceLabels','on','FaceAlpha',0.5)
view(30,30)
title('Bracket with Face Labels')

Figure contains an axes. The axes with title Bracket with Face Labels contains 3 objects of type quiver, patch, line.

figure
pdegplot(model,'FaceLabels','on','FaceAlpha',0.5)
view(-134,-32)
title('Bracket with Face Labels, Rear View')

Figure contains an axes. The axes with title Bracket with Face Labels, Rear View contains 3 objects of type quiver, patch, line.

Установите граничное условие Дирихле со значением 1000 на задней стороне, которая является поверхностью 4. Установите большие поверхности 1 и 7, и также круговую поверхность 11, чтобы иметь Неймановы граничные условия со значением g = -10. Не делайте условий границы множества на других поверхностях. Те поверхности значение по умолчанию к Неймановым граничным условиям со значением g = 0.

applyBoundaryCondition(model,'Face',4,'u',1000);
applyBoundaryCondition(model,'Face',[1,7,11],'g',-10);

Установите c коэффициент к 1, f к 0,1, и a к нелинейному значению '0.1 + 0.001*u.^2'.

c = 1;
f = 0.1;
a = '0.1 + 0.001*u.^2';

Сгенерируйте mesh и решите УЧП. Запустите с исходного предположения u0 = 1000, который совпадает со значением, вы устанавливаете на поверхности 4. Включите Report опция, чтобы наблюдать сходимость во время решения.

generateMesh(model);
u = pdenonlin(model,c,a,f,'U0',1000,'Report','on');
Iteration     Residual     Step size  Jacobian: full
   0          7.2059e-01
   1          1.3755e-01   1.0000000
   2          4.0800e-02   1.0000000
   3          1.1343e-02   1.0000000
   4          2.2740e-03   1.0000000
   5          1.7779e-04   1.0000000
   6          1.4223e-06   1.0000000

Постройте решение на контуре геометрии.

pdeplot3D(model,'ColorMapData',u)

Входные параметры

свернуть все

Модель PDE в виде PDEModel объект.

Пример: model = createpde

Коэффициент УЧП в виде скаляра, матрицы, вектора символов, символьного массива, строкового скаляра, представляет в виде строки вектор или коэффициентную функцию. c представляет коэффициент c в скалярном УЧП

(cu)+au=f

или в системе УЧП

(cu)+au=f

Пример: 'cosh(x+y.^2)'

Типы данных: double | char | string | function_handle
Поддержка комплексного числа: Да

Коэффициент УЧП в виде скаляра, матрицы, вектора символов, символьного массива, строкового скаляра, представляет в виде строки вектор или коэффициентную функцию. a представляет коэффициент a в скалярном УЧП

(cu)+au=f

или в системе УЧП

(cu)+au=f

Пример: 2*eye(3)

Типы данных: double | char | string | function_handle
Поддержка комплексного числа: Да

Коэффициент УЧП в виде скаляра, матрицы, вектора символов, символьного массива, строкового скаляра, представляет в виде строки вектор или коэффициентную функцию. f представляет коэффициент f в скалярном УЧП

(cu)+au=f

или в системе УЧП

(cu)+au=f

Пример: char('sin(x)';'cos(y)';'tan(z)')

Типы данных: double | char | string | function_handle
Поддержка комплексного числа: Да

Граничные условия в виде граничной матрицы или массива данных граничных условий. Передайте массив данных граничных условий как указатель на функцию или как имя файла. Граничная матрица обычно является экспортом из приложения PDE Modeler.

Пример: b = 'circleb1', b = "circleb1", или b = @circleb1

Типы данных: double | char | string | function_handle

Поймайте в сети точки в виде 2 Np матрица точек, где Np число точек в mesh. Для описания (pET) матрицы, смотрите Данные о Mesh, когда [p, e, t] Утраивается.

Как правило, вы используете pE, и t данные, экспортированные из приложения PDE Modeler или сгенерированные initmesh или refinemesh.

Пример: [p,e,t] = initmesh(gd)

Типы данных: double

Поймайте в сети ребра в виде 7- Ne матрица ребер, где Ne количество ребер в mesh. Для описания (pET) матрицы, смотрите Данные о Mesh, когда [p, e, t] Утраивается.

Как правило, вы используете pE, и t данные, экспортированные из приложения PDE Modeler или сгенерированные initmesh или refinemesh.

Пример: [p,e,t] = initmesh(gd)

Типы данных: double

Поймайте в сети треугольники в виде 4- Nt матрица треугольников, где Nt количество треугольников в mesh. Для описания (pET) матрицы, смотрите Данные о Mesh, когда [p, e, t] Утраивается.

Как правило, вы используете pE, и t данные, экспортированные из приложения PDE Modeler или сгенерированные initmesh или refinemesh.

Пример: [p,e,t] = initmesh(gd)

Типы данных: double

Аргументы в виде пар имя-значение

Задайте дополнительные разделенные запятой пары Name,Value аргументы. Name имя аргумента и Value соответствующее значение. Name должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

Пример: 'Jacobian','full'

Приближение якобиана в виде 'full', 'fixed', или 'lumped'.

  • 'full' средняя численная оценка полного якобиана на основе разреженной версии numjac функция. 3-D геометрия использует только 'full', любая другая спецификация дает к ошибке.

  • 'fixed' задает матрицу итерации фиксированной точки, где якобиан аппроксимирован матрицей жесткости. Это - 2D значение по умолчанию геометрии.

  • 'lumped' задает “смешанное” приближение как описано в Нелинейных уравнениях. Это приближение основано на числовом дифференцировании коэффициентов.

Пример: u = pdenonlin(model,c,a,f,'Jacobian','full')

Типы данных: char | string

Начальное предположение решения в виде скаляра, вектора из символов или вектора из чисел. Скаляр задает постоянное начальное условие или для скаляра или для системы УЧП.

Для систем уравнений N и mesh с Np узлы, дайте вектор-столбец с N *Np компоненты. Узлами является любой model.Mesh.Nodes, или p данные из initmesh или meshToPet. Смотрите данные о mesh, когда [p, e, t] утраивается.

Первый Np элементы содержат значения 1 компонента, где значение элемента k соответствует узлу p(k). Следующий Np точки содержат значения 2 компонента и т.д. Может быть удобно сначала представлять начальные условия u0 как Np- N матрица, где первый столбец содержит записи для 1 компонента, второй столбец, содержит записи для 2 компонента и т.д. Итоговым представлением начальных условий является u0(:).

Пример: u = pdenonlin(model,c,a,f,'U0','x.^2-y.^2')

Типы данных: double | char | string
Поддержка комплексного числа: Да

Остаточный размер при завершении в виде положительной скалярной величины. pdenonlin выполняет итерации, пока остаточный размер не меньше 'Tol'.

Пример: u = pdenonlin(model,c,a,f,'Tol',1e-6)

Типы данных: double

Максимальное количество итераций Ньютона Гаусса в виде положительного целого числа.

Пример: u = pdenonlin(model,c,a,f,'MaxIter',12)

Типы данных: double

Минимальное затухание поискового направления в виде положительной скалярной величины.

Пример: u = pdenonlin(model,c,a,f,'MinStep',1e-3)

Типы данных: double

Распечатайте информацию о сходимости в виде 'off' или 'on'.

Пример: u = pdenonlin(model,c,a,f,'Report','on')

Типы данных: char | string

Норма невязки в виде p значение для нормы LP, или как 'energy'P может быть любое положительное действительное значение, Inf, или -Inf. p норма векторного v sum(abs(v)^p)^(1/p). Смотрите norm.

Пример: u = pdenonlin(model,c,a,f,'Norm',2)

Типы данных: double | char | string

Выходные аргументы

свернуть все

Решение для УЧП, возвращенное как вектор.

  • Если УЧП является скаляром, означая только одно уравнение, то u вектор-столбец, представляющий решение u в каждом узле в mesh. u(i) решение в iстолбец th model.Mesh.Nodes или iстолбец th p.

  • Если УЧП является системой N> 1 уравнение, то u вектор-столбец с N *Np элементы, где Np количество узлов в mesh. Первый Np элементы u представляйте решение уравнения 1, затем следующий Np элементы представляют решение уравнения 2 и т.д.

Чтобы получить решение в произвольной точке в геометрии, используйте pdeInterpolant.

Чтобы построить решение, использовать pdeplot для 2D геометрии, или см. 3-D Графики Решения и Градиента с MATLAB® Functions.

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

Советы

  • Если итерация Ньютона не сходится, pdenonlin отображает сообщение об ошибке Too many iterations или Stepsize too small.

  • Если исходное предположение производит матрицы, содержащие NaN или Inf элементы, pdenonlin отображает сообщение об ошибке Unsuitable initial guess U0 (default: U0 = 0).

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

Алгоритмы

свернуть все

Нелинейные уравнения

Основная идея состоит в том, чтобы использовать итерации Ньютона Гаусса, чтобы решить нелинейные уравнения. Скажите, что вы пытаетесь решить уравнение

r (u) = – ∇ · (c (u) ∇u) + a (u) u - f (u) = 0.

В установке FEM вы решаете слабую форму r (u) = 0. Установите, как обычно,

u(x)=Ujϕj

где x представляет 2D или 3-D точку. Затем умножьте уравнение на произвольную тестовую функцию ϕi, объединяйтесь на области Ω и используйте формулу Грина и граничные условия, чтобы получить

0=ρ(U)=j(Ω((c(x,U)ϕj(x))ϕj(x)+a(x,U)ϕj(x)ϕi(x))dx+Ωq(x,U)ϕj(x)ϕi(x)ds)UjΩf(x,U)ϕi(x)dxΩg(x,U)ϕi(x)ds

который должен содержать для всех индексов i.

Вектор невязок ρ (U) может быть легко вычислен как

ρ (U) = (K + M + Q) U – (F + G)

где матрицы K, M, Q и векторы F и G производятся путем сборки проблемы

– ∇ · (c (U) ∇u) + a (U) u = f (U).

Примите, что у вас есть предположение U (n) решения. Если U (n) достаточно близок к точному решению, улучшенное приближение, U (n +1) получен путем решения линеаризовавшей задачи

ρ(U(n))U(U(n+1)U(n))=αρ(U(n))

где α положительное число. (Не необходимо, чтобы ρ (U) = 0 имел решение, даже если ρ (u) = 0 имеет.) В этом случае итерация Ньютона Гаусса имеет тенденцию быть минимизатором невязки, т.е. решением minU ρ(U).

Это известно это за достаточно маленький α

ρ(U(n+1))<ρ(U(n))

и

pn=(ρ(U(n))U)1ρ(U(n))

называется направлением спуска для ρ(U), где 2-норма L. Итерация

U(n+1)=U(n)+αpn,

где α ≤ 1 выбран как можно больше таким образом, что шаг имеет разумный спуск.

Метод Ньютона Гаусса локален, и сходимость гарантируют только, когда U (0) достаточно близок к решению. В общем случае первое предположение может найтись вне области сходимости. Чтобы улучшить сходимость от плохих исходных предположений, стратегия затухания реализована для выбора α, поиск линии Армиджо-Голдстайна. Это выбирает самый большой коэффициент демпфирования α из последовательности 1, 1/2, 1/4... таким образом, что следующее неравенство содержит:

ρ(U(n))ρ(U(n))+αpnα2ρ(U(n))

который гарантирует сокращение нормы невязки по крайней мере 1 – α/2. Каждый шаг алгоритма поиска линии требует оценки невязки ρ(U(n)+αpn).

Важный момент этой стратегии - это, когда U (n) приближается к решению, затем α→1 и таким образом повышения ставки сходимости. Если существует решение ρ (U) = 0, схема в конечном счете восстанавливает квадратичный уровень сходимости стандарта итерация Ньютона.

Тесно связанный с предыдущей проблемой выбор исходного предположения U (0). По умолчанию решатель устанавливает U (0) и затем собирает матрицы FEM K и F и вычисляет

U (1) = K –1F

Ослабленная итерация Ньютона Гаусса затем начата с U (1), который должен быть лучшим предположением, чем U (0). Если граничные условия не зависят от решения u, то U (1) удовлетворяет им, даже если U (0) не делает. Кроме того, если уравнение линейно, то U (1) является точным решением FEM, и решатель не вводит цикл Ньютона Гаусса.

Существуют ситуации, где U (0) = 0 не имеет никакого смысла, или сходимость невозможна.

В некоторых ситуациях у вас может уже быть хорошее приближение, и нелинейный решатель может быть начат с него, избежав медленного режима сходимости. Эта идея используется в адаптивном генераторе mesh. Это вычисляет решение U˜ на mesh, оценивает ошибку и может совершенствовать определенные треугольники. interpolant U˜ очень хорошее начальное приближение для решения на усовершенствованной mesh.

В целом точный якобиан

Jn=ρ(U(n))U

не доступно. Приближение Jn конечными разностями следующим образом является дорогим, но выполнимым. i th столбец Jn может быть аппроксимирован

ρ(U(n)+εϕi)ρ(U(n))ε

который подразумевает сборку матриц FEM для треугольников, содержащих узел решетки i. Очень простое приближение к Jn, который дает итерацию фиксированной точки, также возможно следующие. По существу, для данного U (n), вычислите матрицы FEM K и F и установите

U (n +1) = K –1F .

Это эквивалентно аппроксимации якобиана с матрицей жесткости. Действительно, начиная с ρ (U (n)) = KU (n)F, помещая Jn = выражения K

U(n+1)=U(n)Jn1ρ(U(n))=U(n)K1(KU(n)F)=K1F

Во многих случаях уровень сходимости является медленным, но стоимость каждой итерации является дешевой.

Partial Differential Equation Toolbox™ нелинейный решатель также предусматривает компромисс между этими двумя экстремальными значениями. Чтобы вычислить производную отображения UKU, продолжите можно следующим образом. Термин a был не использован для ясности, но появляется снова в конечном результате.

(KU)iUj=limε01εl(Ωc(U+εϕj)ϕlϕidx(Ul+εδl,j)Ωc(U)ϕlϕidxUl)=Ωc(U)ϕjϕidx+lΩϕjcuϕlϕidxUl

Первый интегральный срок является не чем иным как Ki,j.

Второй срок “смешан”, т.е. заменен диагональной матрицей, которая содержит суммы строки. С тех пор Σjϕj = 1, второй срок аппроксимирован

δi,jlΩcuϕlϕidxUl

который является i th компонент K (c') U, где K (c') является матрицей жесткости, сопоставленной с коэффициентом ∂c / ∂ u, а не c. То же обоснование может быть применено к производной отображения UMU. Производная отображения U → –F точно

Ωfuϕiϕjdx

который является большой матрицей, сопоставленной с коэффициентом ∂f / ∂ u. Таким образом якобиан остаточного ρ (U) аппроксимирован

J=K(c)+M(af)+diag((K(c)+M(a))U)

где дифференцирование относительно u, K и M определяют жесткость и большие матрицы, и их индексы определяют коэффициенты, относительно которых они собраны. В каждой итерации Ньютона Гаусса нелинейный решатель собирает матрицы, соответствующие уравнениям

(cu)+(af')u=0(c'u)+a'u=0

и затем производит аппроксимированный якобиан. Дифференцирования коэффициентов сделаны численно.

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

K˜U˜=[KHH0][Uμ]=[FR]=F˜

где коэффициенты K˜ и F˜ может зависеть от решения U˜. “Смешанный” подход аппроксимирует производное отображение невязки

[JHH0]

Нелинейность граничных условий и зависимости коэффициентов на производных U˜ правильно не линеаризуются этой схемой. Когда такая нелинейность сильна, схема уменьшает до итерации фиксировать-точки и может сходиться медленно или нисколько. Когда граничные условия линейны, они не влияют на свойства сходимости схем итерации. В Неймановом случае они невидимы (H является пустой матрицей), и в случае Дирихле они просто утверждают, что невязка является нулем на соответствующих граничных точках.

Смотрите также

Представлено до R2006a
Для просмотра документации необходимо авторизоваться на сайте