Ctrl+K

Пример. Нелинейное уравнение Пуассона

Пример. Нелинейное уравнение Пуассона#

Или как работать с нелинейностями?

\[\begin{split} \begin{cases}-\nabla \cdot(q(u) \nabla u)=0 & \text { in } \Omega=[0,1], \\ u=u_0 & \text { on } \partial \Omega\end{cases} \end{split}\]

Есть два простых метода решения этой задачи. Оба основаны на последовательном приближении к ответу.

1. Метод итераций (метод Пикарда)#

Идея: по текущему прилижению определить следующее

\[ -\nabla \cdot\left(q^k\left(u^k\right) \nabla u^{k+1}\right)=0 \quad \text { in } \Omega=[0,1], \]

В слабой форме получим:

\[ \int_{\Omega} q\left(u^k\right) \nabla u^{k+1} \cdot \nabla v \mathrm{~d} x=0 \]

Продолжаем итерацию, пока решение на новом шаге будет слабо отличатся от текущего решения.

2. Метод Ньютона#

Идея: по текущему приближению определить возмущение, которое улучшит приближение.

\[\begin{split} \begin{aligned} & u^{k+1}=u^k+\delta u \\ & q\left(u^{k+1}\right)=q\left(u^k\right)+q^{\prime}\left(u^k\right) \delta u+\mathcal{O}\left((\delta u)^2\right) \approx q\left(u^k\right)+q^{\prime}\left(u^k\right) \delta u \end{aligned} \end{split}\]

Подставим это в исходное уравнение, отбрасывая члены второго порядка.

\[\begin{split} \begin{gathered} \nabla \cdot\left(q\left(u^k\right) \nabla u^k\right)+\nabla \cdot\left(q\left(u^k\right) \nabla \delta u\right)+\nabla \cdot\left(q^{\prime}\left(u^k\right) \delta u \nabla u^k\right)=0 \\ \nabla \cdot\left(q\left(u^k\right) \nabla \delta u\right)+\nabla \cdot\left(q^{\prime}\left(u^k\right) \delta u \nabla u^k\right)=-\nabla \cdot\left(q\left(u^k\right) \nabla u^k\right) \end{gathered} \end{split}\]

Домножая на тестовую функцию и проводя интегрирование по частям, приходим к слабой форме.

\[ \int_{\Omega}\left(q\left(u^k\right) \nabla \delta u \cdot \nabla v+q^{\prime}\left(u^k\right) \delta u \nabla u^k \cdot \nabla v\right) \mathrm{d} x=-\int_{\Omega} q\left(u^k\right) \nabla u^k \cdot \nabla v \mathrm{~d} x \]

Видно, что уравнение имеет вид:

Билинейная форма \((\delta u, v)=\) Линейная форма \((v)\).

Важные детали:

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

  2. Возмущение должно зануляться на границе

mesh = fe.UnitIntervalMesh(40)
V = fe.FunctionSpace(mesh, "Lagrange", 1)


tol = 1e-14

def left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < tol


def right_boundary(x, on_boundary):
    return on_boundary and abs(x[0]-1) < tol


Gamma_0 = fe.DirichletBC(V, fe.Constant(0.0), left_boundary)
Gamma_1 = fe.DirichletBC(V, fe.Constant(1.0), right_boundary)


bcs = [Gamma_0, Gamma_1]
# Найдем начальное приближение (q(u)=1, m=0)
u = fe.TrialFunction(V)
v = fe.TestFunction(V)
a = fe.inner(fe.nabla_grad(u), fe.nabla_grad(v))*fe.dx
f = fe.Constant(0.0)
L = f*v*fe.dx
A, b = fe.assemble_system(a, L, bcs)
u_k = fe.Function(V)
U_k = u_k.vector()
fe.solve(A, U_k, b)

1

fe.plot(u_k)
plt.title('Inital guess')

Text(0.5, 1.0, 'Inital guess')
../../_images/30922650e0f951cb7ba79342f555b5a11301ffd7ab17c3fa89ba6091d91ca2a6.png 
#добавим граничные условия для возмущения du
Gamma_0_du = fe.DirichletBC(V, fe.Constant(0), left_boundary)
Gamma_1_du = fe.DirichletBC(V, fe.Constant(0), right_boundary)
bcs_du = [Gamma_0_du, Gamma_1_du]

m = 8

def q(u):
    return (1+u)**m

#производная от q
def Dq(u):
    return m*(1+u)**(m-1)


du = fe.TrialFunction(V) # u = u_k + omega*du
a = fe.inner(q(u_k)*fe.nabla_grad(du), fe.nabla_grad(v))*fe.dx + \
    fe.inner(Dq(u_k)*du*fe.nabla_grad(u_k), fe.nabla_grad(v))*fe.dx
L = -fe.inner(q(u_k)*fe.nabla_grad(u_k), fe.nabla_grad(v))*fe.dx

du = fe.Function(V)
u = fe.Function(V) # u = u_k + omega*du
omega = 1.0 # relaxation parameter
eps = 1.0
tol = 1.0E-5
iter = 0
maxiter = 25
while eps > tol and iter < maxiter:
    iter += 1
    A, b = fe.assemble_system(a, L, bcs_du)
    fe.solve(A, du.vector(), b)
    eps = np.linalg.norm(du.vector(), ord=np.Inf)
    print("Norm:", eps)
    u.vector()[:] = u_k.vector() + omega*du.vector()
    u_k.assign(u)

Norm: 2.6926580821284225
Norm: 0.42678052234021463
Norm: 0.3790764929916502
Norm: 0.3362295391794988
Norm: 0.29701382866741033
Norm: 0.25931623334952986
Norm: 0.2189422903508294
Norm: 0.1703542478525539
Norm: 0.10525359212091752
Norm: 0.03666912499970038
Norm: 0.0037568546653569818
Norm: 3.590420381359272e-05
Norm: 3.2438027096345295e-09

fe.plot(u_k)

[<matplotlib.lines.Line2D at 0x7c4deb1707c0>]
../../_images/42971fe4d277d7d3e6fa827b57f60f194d4e6271359ce7f9d8afee41b1b69d8e.png
Содержание