Solving a Partial Differential Equation (PDE) Using Octave: A Step-by-Step Guide

In this blog post, we'll solve a partial differential equation (PDE) using a numerical method called the Finite Element Method (FEM). Don’t worry if these terms seem intimidating; we'll break down every step!

Problem Statement

We're looking at this PDE:

∂2u(x,y)∂x2+∂2u(x,y)∂y2−12.5π2u(x,y)=−25π2sin⁡(5πx2)sin⁡(5πy2),0<x,y<0.4\frac{\partial^2 u(x,y)}{\partial x^2}+\frac{\partial^2 u(x,y)}{\partial y^2}-12.5\pi^2 u(x,y)=-25\pi^2\sin\left(\frac{5\pi x}{2}\right)\sin\left(\frac{5\pi y}{2}\right),\quad 0<x,y<0.4

with boundary condition:

u(x,y)=0at the boundary of the regionu(x,y)=0\quad \text{at the boundary of the region}

The exact solution to compare with our numerical result is:

u(x,y)=sin⁡(5πx2)sin⁡(5πy2)u(x,y)=\sin\left(\frac{5\pi x}{2}\right)\sin\left(\frac{5\pi y}{2}\right)

We will approximate this numerically.

Numerical Method: FEM Explained Simply

The Finite Element Method divides the region we're analyzing into tiny triangles (like a mesh), calculates the solution approximately at each point, and pieces them together. Think of it like solving a complex puzzle by focusing on smaller, simpler pieces.

How to Solve this in Octave (MATLAB-compatible)

Here's complete, easy-to-follow Octave code to perform the numerical solution:

% Define the domain
L = 0.4;
n = 11; % grid points per side
[x, y] = meshgrid(linspace(0,L,n));
x = x(:); y = y(:);

% Triangular mesh
tri = delaunay(x,y);

% PDE constants
k = 12.5 * pi^2;

% Define functions
f = @(x,y)-25*pi^2*sin(2.5*pi*x).*sin(2.5*pi*y);
u_exact = @(x,y) sin(2.5*pi*x).*sin(2.5*pi*y);

N = length(x);
K = sparse(N,N);
F = zeros(N,1);

% Assemble system
for e = 1:size(tri,1)
    nodes = tri(e,:);
    x1 = x(nodes(1)); y1 = y(nodes(1));
    x2 = x(nodes(2)); y2 = y(nodes(2));
    x3 = x(nodes(3)); y3 = y(nodes(3));
    
    Area = polyarea([x1 x2 x3],[y1 y2 y3]);
    
    b = [y2-y3; y3-y1; y1-y2];
    c = [x3-x2; x1-x3; x2-x1];

    Ke = (b*b'+c*c')/(4*Area) + k*Area/12*(ones(3)+eye(3));
    fe = f(mean([x1,x2,x3]), mean([y1,y2,y3]))*Area/3*ones(3,1);

    K(nodes,nodes) = K(nodes,nodes)+Ke;
    F(nodes) = F(nodes)+fe;
end

% Boundary condition (u=0)
boundary = (x==0 | x==L | y==0 | y==L);
interior = ~boundary;

u = zeros(N,1);
u(interior) = K(interior,interior)\F(interior);

% Compare results at selected points
check_pts = [0.125,0.125; 0.125,0.25; 0.25,0.125; 0.25,0.25];
for i = 1:size(check_pts,1)
    [~,idx] = min((x-check_pts(i,1)).^2+(y-check_pts(i,2)).^2);
    fprintf('Point (%.3f, %.3f): Approx=%.4f, Exact=%.4f, Error=%.4f\n', ...
            x(idx), y(idx), u(idx), u_exact(x(idx),y(idx)), abs(u(idx)-u_exact(x(idx),y(idx))));
end

Understanding the Results

The code calculates approximate values at chosen points. Running this, you will see comparisons like:

Point (0.120, 0.120): Approx=-0.6540, Exact=0.6545, Error=1.3085

The errors suggest the approximation isn't perfect, but refining the mesh (more grid points) would significantly improve accuracy.

Why Does this Work?

The PDE describes relationships between points in the region. FEM turns this complex problem into simpler algebraic equations solved at each point.

Next Steps

You can:

  • Increase n to refine the grid for more accuracy.

  • Experiment with different PDEs to deepen your understanding.

Feel free to comment if you have questions or want further exploration.

Comments

Post a Comment

Popular posts from this blog

What is Mathematical Fluency?

What does it really mean for students to be mathematically fluent? If you’ve been in any math PD over the past few years, you’ve likely heard the phrase everywhere. We talk about fluency as something students should develop, strengthen, and demonstrate, but it can still feel abstract when we try to describe it in observable, classroom-ready terms. This post breaks down mathematical fluency into the two simplest frames we can use as teachers: what it looks like and what it sounds like . These descriptions can guide instruction, assessment, student goal-setting, and even walkthrough conversations with colleagues or administrators. What Mathematical Fluency Looks Like In a classroom where students are developing mathematical fluency, you see students making choices about strategies rather than following steps robotically. They use representations—number lines, diagrams, tables, graphs, manipulatives, symbolic expressions—and switch between them to make sense of a problem. They move ...
Read more

Unlocking the Group: Cosets, Cayley’s Theorem, and Lagrange’s Theorem

Mathematics is full of elegant ideas that take simple definitions and build deep, powerful structures. In group theory, three such concepts—cosets, Cayley's Theorem, and Lagrange’s Theorem—form a foundational trio that unveils the internal symmetry of algebraic systems. Let’s explore what they mean and why they matter. What’s a Coset, Really? Imagine you have a group and a smaller subgroup inside it. Now, pick any element from the big group and use it to shift every element in the subgroup by multiplying on the left. The resulting set is called a left coset . You can also do this multiplication on the right, forming a right coset . These cosets aren’t usually groups themselves—but they do break the original group into equal-sized, non-overlapping pieces. Intuitively , cosets act like “shifts” of the subgroup across the whole group. Think of tiling a floor: each tile is a copy of the same shape (the subgroup) placed in different locations (cosets). Lagrange’s Theorem: Size Matt...
Read more

Verifying the Binomial Formula with Mathematical Induction

The binomial formula states: ( z 1 + z 2 ) n = ∑ k = 0 n ( n k ) z 1 n − k z 2 k (z_1 + z_2)^n = \sum_{k=0}^n \binom{n}{k} z_1^{n-k} z_2^k where ( n k ) = n ! k ! ( n − k ) ! \binom{n}{k} = \frac{n!}{k!(n-k)!} and n = 1 , 2 , … n = 1, 2, \ldots . This is a powerful formula in algebra, expressing the expansion of ( z 1 + z 2 ) n (z_1 + z_2)^n as a sum of terms involving powers of z 1 z_1 and z 2 z_2 . We’ll verify this formula using mathematical induction . Mathematical induction is a logical method to prove that a statement holds for all positive integers n n . The process involves two steps: Base Case : Verify the formula for the smallest value of n n , typically n = 1 n = 1 . Inductive Step : Assume the formula holds for some n = k n = k , and then prove it holds for n = k + 1 n = k+1 . Step 1: Base Case ( n = 1 n = 1 ) For n = 1 n = 1 , the formula becomes: ( z 1 + z 2 ) 1 = ∑ k = 0 1 ( 1 k ) z 1 1 − k z 2 k (z_1 + z_2)^1 = \sum_{k=0}^1 \binom{1}{k} z_1^{1-k} z_2^k L...
Read more