Finite Element Method for Gravitational Potential Equation

Project Overview

This project implements the Finite Element Method (FEM) to solve the gravitational potential differential equation:

d²Φ/dx² = 4πGρ(x)

with boundary conditions:

• Φ(0) = 5
• Φ(3) = 4

and density function:

• ρ(x) = 0 for x ∈ [0,1]
• ρ(x) = 1 for x ∈ (1,2]
• ρ(x) = 0 for x ∈ (2,3]

Project Structure

Files:

• wyprowadzenie.pdf - Mathematical derivation of the variational formulation and FEM implementation
• fem.m - MATLAB implementation of the FEM solver

Mathematical Background

Weak Formulation

The strong form is converted to weak formulation using:

• Test functions v ∈ V that vanish at boundaries
• Integration by parts to reduce derivative order
• Construction of bilinear form B(Φ, v) and linear functional L(v)

Shift Function

Due to non-homogeneous Dirichlet boundary conditions, a shift function is introduced: Φ(x) = w(x) + Φ̃(x) where Φ̃(x) = -1/3*x + 5 satisfies the boundary conditions.

Finite Element Discretization

• Domain Ω = [0,3] divided into N equal elements
• Linear basis functions e_i(x) with compact support
• Piecewise linear approximation of the solution

Usage

Parameters:

• N = 1000 - Number of elements (adjustable)
• G = 6.6743e-11 - Gravitational constant (adjustable)

Running the Code:

% Simply run the fem.m file in MATLAB
% The code will:
% 1. Set up the finite element discretization
% 2. Assemble the stiffness matrix and load vector
% 3. Solve the linear system
% 4. Plot the solution Φ(x)

Output:

• Graphical plot of the gravitational potential Φ(x) over domain [0,3]
• Numerical solution values stored in variable y

Requirements

• MATLAB
• No additional toolboxes required

Features

• Adjustable number of elements N
• Configurable gravitational constant G
• Numerical integration using Gauss-Legendre quadrature
• Automatic handling of non-homogeneous boundary conditions
• Visualization of the solution

Numerical Method

• Method: Finite Element Method (FEM)
• Element Type: Linear Lagrangian elements
• Integration: 2-point Gauss-Legendre quadrature
• Matrix Solver: MATLAB's backslash operator (direct solver)

This implementation demonstrates the complete FEM workflow from mathematical formulation to numerical solution for a second-order boundary value problem with piecewise constant coefficients.