Potencjał grawitacyjny MES - Szymon Iwaniuk

Mix.install([
  {:nx, "~> 0.7"},
  {:vega_lite, "~> 0.1.8"},
  {:kino_vega_lite, "~> 0.1.10"}
])

1. Treść problemu

$$ \frac{d^{2}\Phi}{dx^{2}} = 4\pi G\rho(x) $$ $$ \Phi(0) = 5 $$ $$ \Phi(3) = 7 $$

$$ \rho(x) = \begin{cases} 0 & \text{dla } x \in [0, 1] \ 1 & \text{dla } x \in [1, 2] \ 0 & \text{dla } x \in (2, 3] \end{cases} $$

$$ G = 1 $$

Szukaną funkcją jest $\Phi(x)$

2. Sformułowanie wariacyjne problemu

$$ \Phi^{\prime\prime}=4\pi G\rho $$ Niech $v:[0,3]\rightarrow\mathbb{R}, v\in V$ gdzie $V$ jest przestrzenią funkcji zerujących się na brzegach. Obustronnie mnożymy przez $v$ powyższe równanie.

$$ \Phi^{\prime\prime}v=4\pi G\rho v $$

Następnie całkujemy je na $[0, 3]$

$$ \int{0}^{3}\Phi^{\prime\prime}v~dx=4\pi G\int{0}^{3}\rho v~dx $$

Stosujemy całkowanie przez części aby uzyskać

$$ [\Phi^{\prime}v]{0}^{3}-\int{0}^{3}\Phi^{\prime}v^{\prime}dx=4\pi G\int_{0}^{3}\rho v~dx $$

Ponieważ $$ v\in V \Rightarrow v(3)=v(0)=0 \Rightarrow [\Phi^{\prime}v]_{0}^{3}=0 $$

Mamy więc

$$ -\int{0}^{3}\Phi^{\prime}v^{\prime}dx=4\pi G\int{0}^{3}\rho v~dx $$

Definiujemy:

$$ B(\Phi,v)=-\int_{0}^{3}\Phi^{\prime}v^{\prime}dx $$

$$ L(v)=4\pi G\int_{1}^{2}v~dx $$

Nasze początkowe równanie ma zatem postać

$$ B(\Phi,v)=L(v) $$

3. Warunki brzegowe

Przez niezerowe warunki Dirichleta musimy skonstruować przesunięcie: $$ \Phi = w + \overline{\Phi} \quad \text{gdzie } w \in V $$

$$ \overline{\Phi} = Ax + B, \quad \overline{\Phi}(0)=5, \quad \overline{\Phi}(3)=7 $$

$$ \begin{cases} 0 \cdot A + B = 5 \ 3 \cdot A + B = 7 \end{cases} \Rightarrow \begin{cases} B = 5 \ A = \frac{2}{3} \end{cases} $$

$$ \overline{\Phi}(x) = 5 + \frac{2}{3}x $$

Tak skonstruowana funkcja $\overline{\Phi}$ spełnia warunki brzegowe, wracając do równania, możemy zapisać je w postaci: $$ B(w+\overline{\Phi}, v) = L(v) $$

B jest liniowe względem pierwszego argumentu, zatem: $$ B(w, v) + B(\overline{\Phi}, v) = L(v) $$

$$ B(w, v) = L(v) - B(\overline{\Phi}, v) $$

Przyjmując nowe oznaczenie $\overline{L}(v) = L(v) - B(\overline{\Phi}, v)$

$$ B(w, v) = \overline{L}(v) $$

4. Konstrukcja podprzestrzeni N elementów skończonych

Konstruujemy przestrzeń $V_h \subset V$, dzieląc przedział $[0,3]$ na $N$ równych części, każda o długości $h = \frac{3}{N}.$ Dla takiego podziału punkty brzegowe to: $$ x_0 = 0,\quad x_i = ih,\quad x_N = 3. $$

$$ x{i-1} = x_i - h,\quad x{i+1} = x_i + h $$

Wybieramy funkcje bazowe w postaci:

$$ ei(x) = \begin{cases} \frac{x - x{i-1}}{h} & \text{dla}\ x \in [x{i-1}, x_i], \ \frac{x{i+1} - x}{h} & \text{dla}\ x \in (xi, x{i+1}], \text{gdzie}\ i\in{1,2,\dots,N-1} \ 0 & \text{w pozostałych przypadkach} \end{cases} $$

Równowanznie:

$$ ei(x) = \begin{cases} \frac{x}{h} - i + 1 & \text{dla}\ x \in [x{i-1}, xi], \ i + 1 - \frac{x}{h} & \text{dla}\ x \in (x_i, x{i+1}], \ 0 & \text{w pozostałych przypadkach} \end{cases} $$

Pochodna funkcji bazowej:

$$ \frac{d(ei)}{dx} = \begin{cases} \frac{1}{h} & \text{dla}\ x \in [x{i-1}, xi], \ -\frac{1}{h} & \text{dla}\ x \in (x_i, x{i+1}], \ 0 & \text{w pozostałych przypadkach} \end{cases} $$

Mozemy wiec zdefiniowac $Vh = \langle e_1, e_2, \dots, e{N-1} \rangle$

5. Problem przybliżony

Z rozwazanego powyzej problemu mozemy przejsc na problem przyblizony. Znaleźć $w{h}\in V{h}$, które spełnia $$ B(w{h},v{h})=\overline{L}(v{h}) \text{ dla każdego } v{h}\in V_{h} $$

$$ w \approx w{h} = \sum{i=1}^{N-1}w{i}e{i} $$

Układ $N-1$ równań przyjmuje postać: $$ B(\sum_{i=1}^{N-1} w_ie_i,v_j) = \overline{L}(v_j) \quad \text{dla } j \in {0,1,\dots,N} $$

Przyjmujemy $v_j = e_j$

$$ B(\sum_{i=1}^{N-1} w_ie_i,e_j) = \overline{L}(e_j) $$

B jest liniowe wzgledem pierwszego argumentu, zatem:

$$ \sum_{i=1}^{N-1} w_iB(e_i,e_j) = \overline{L}(e_j) $$

Uklad równań w postaci macierzowej:

$$ \begin{pmatrix} B(e1,e_1) & B(e_2,e_1) & \cdots & B(e{N-1},e1) \ B(e_1,e_2) & B(e_2,e_2) & \cdots & B(e{N-1},e2) \ \vdots & \vdots & \ddots & \vdots \ B(e_1,e{N-1}) & B(e2,e{N-1}) & \cdots & B(e{N-1},e{N-1}) \end{pmatrix} * \begin{pmatrix} w1 \ w_2 \ \vdots \ w{N-1} \end{pmatrix}

\begin{pmatrix} \overline{L}(e1) \ \overline{L}(e_2) \ \vdots \ \overline{L}(e{N-1}) \end{pmatrix} $$

Musimy zatem policzyć

$$ B(e{i},e{j})=-\int{0}^{3}e{i}^{\prime}e_{j}^{\prime}dx $$

$$ \overline{L}(e{j})=L(e{j})-B(\overline{\Phi},e_{j}) $$

$$ L(e{j})=4\pi G\int{1}^{2}e_{j}dx $$

Łatwo zauwazyć, ze $B(e{i},e{j})=0$ gdy $j < i - 1$ lub $j > i + 1$. W przeciwnym wypadku $$ B(ei,e{i+1}) = B(e{i+1},e_i) = -\int{xi}^{x{i+1}} ei’\, e{i+1}’\, dx. $$

$$ B(ei,e_i) = -\int{x{i-1}}^{x{i+1}} e_i’\, e_i’\, dx, \quad \text{dla } i \in {1,2,\dots,N-1}. $$

Zauwazamy rowniez, ze $L(ej) = 0$ gdy $hj < 1$ lub $hj > 2$. W przeciwnym wypadku: $$ L(e_j) = 4\pi G\int{x{j-1}}^{x{j+1}} e_j\ dx, \quad \text{dla } j \in {1,2,\dots,N-1}. $$

Ostatecznie, pamietając ze ${\Phi}’ = \frac{2}{3}$ $$ B(\overline{\Phi},ej) = -\frac{2}{3} \int{x{j-1}}^{x{j+1}} e_j’\ dx \quad \text{dla } j \in {1,2,\dots,N-1}. $$

6. Kwadratura Gaussa

Do przybliżenia wartości całek stosujemy kwadraturę Gaussa-Legendre’a, zamieniając przedziały całkowania na $[-1, 1]$: $$ \int{a}^{b} f(x) dx = \int{-1}^{1}f(\frac{b-a}{2}t + \frac{a + b}{2})\frac{b-a}{2}dt $$

Wykorzystując, ze

$$ \int{-1}^{1}f(x)dx \approx\sum{i=1}^{N}w_if(x_i) $$

Dla przedziału $[a, b]$ otrzymujemy:

$$ \int{a}^{b}f(x)dx\approx\frac{b-a}{2}\sum{i=1}^{n}w{i}f(\frac{b-a}{2}t{i}+\frac{a+b}{2}) $$

Przyjmujemy $n=2 \Rightarrow w{i}=1$ $ \wedge$ $t{i}=\pm\frac{1}{\sqrt{3}}$ $$ \int_{a}^{b}f(x)dx\approx \frac{b-a}{2}\left [f\left(\frac{b-a}{2}\frac{1}{\sqrt{3}}+\frac{a+b}{2}\right)+ f\left(\frac{a-b}{2}\frac{1}{\sqrt{3}}+\frac{a+b}{2}\right)\right] $$

7. Rozwiązanie

Ostatecznie rozwiązanie otrzymamy z pomocą programu w języku Elixir, który obliczy wartości $w0,w_1,\dots,w_N$. Dzięki tym wartościom uzyskamy: $$ {\Phi} = w + \overline{\Phi} \approx \sum{i=0}^{N}w_i e_i + 5 + \frac{2x}{3} $$

_metadata = 
"""
Created on Wed Dec 26 12:04:45 2025 

@author: Szymon Iwaniuk
"""

# 1. PARAMETERS OF THE MODEL
# choose different parameters based on your needs
g = 1.0           # Gravitational const
l = 3.0           # Total length of the domain
phi_0 = 5.0       # Potential at x = 0
phi_l = 7.0       # Potential at x = 3

# 2. DISCRETIZATION
n = 50            # Number of finite elements 
h = l / n         # Width of each element 

# 3. MODEL
# Denisty function
rho = fn x ->
  cond do
    0.0 <= x and x <= 1.0 -> 0.0
    1.0 < x and x <= 2.0 -> 1.0
    2.0 < x and x <= 3.0 -> 0.0
    true -> raise "#{x} out of bounds"
  end 
end

# Basis function
e_i = fn x, i, h ->
  xi_minus_1 = (i - 1) * h
  xi = i * h
  xi_plus_1 = (i + 1) * h

  cond do
    x >= xi_minus_1 and x <= xi -> 
      x / h - i + 1
    x > xi and x <= xi_plus_1 -> 
      i + 1 - x / h
    true -> 
      0.0
  end
end

# Derivative of basis function 
de_i = fn x, i, h ->
  xi_minus_1 = (i - 1) * h
  xi = i * h
  xi_plus_1 = (i + 1) * h

  cond do
    x >= xi_minus_1 and x <= xi -> 
      1 / h
    x > xi and x <= xi_plus_1 -> 
      - 1 / h
    true -> 
      0.0
  end
end

# 5. INITIALIZATION OF THE GRID OF VALUES
# Const before integral
coeff = 4 * :math.pi() * g # 4 * pi * G

# Create matrix B
matrix_B_list = for i <- 0..n do
  for j <- 0..n do
    cond do
      i == j -> -2.0 / h
      abs(i - j) == 1 -> 1.0 / h
      true -> 0.0
    end
  end
end

# Create vector L
vector_L_list = for i <- 0..n do
  xi = i * h
  rho.(xi) * coeff * h
end

# 6. BOUNDARIES
final_matrix = Enum.with_index(matrix_B_list) |> Enum.map(fn {row, i} ->
  cond do
    i == 0 -> [1.0 | List.duplicate(0.0, n)] # ϕ(0) = phi_0
    i == n -> List.duplicate(0.0, n) ++ [1.0] # ϕ(l) = phi_l
    true -> row
  end
end)

final_vector = Enum.with_index(vector_L_list) |> Enum.map(fn {val, i} ->
  cond do
    i == 0 -> phi_0
    i == n -> phi_l
    true -> val
  end
end)

# Conversion to Nx tensors
b_final = Nx.tensor(final_matrix)
l_final = Nx.tensor(final_vector)

# 7. SOLVING THE SYSTEM
phi_results = Nx.LinAlg.solve(b_final, l_final)

# 8. PLOTTING THE SURFACE
phi_list = Nx.to_flat_list(phi_results)
plot_data = Enum.with_index(phi_list) 
            |> Enum.map(fn {p, i} -> %{"x" => i * h, "ϕ" => p} end)

alias VegaLite, as: Vl
Vl.new(width: 800, height: 600, title: "y = ϕ(x)")
|> Vl.data_from_values(plot_data)
|> Vl.mark(:line, point: true)
|> Vl.encode_field(:x, "x", type: :quantitative,
  axis: [title_font_size: 20])
|> Vl.encode_field(:y, "ϕ", type: :quantitative,
  axis: [title_font_size: 20])