ODE and PDE implementations in MATLAB

mandatory_2.m

%#ok<*NOPTS>
format shortG
clear; close all; clc;

% two coupled equations with integral terms:

% u(x) = (e1 * sigma * T1^4) + (1 - e1) * ∫[-0.5w;0.5w] F(x, y, d) * v(y) dy
% v(y) = (e2 * sigma * T2^4) + (1 - e2) * ∫[-0.5w;0.5w] F(x, y, d) * u(x) dx

% constants
global T1 T2 e1 e2 sigma d w
[T1, T2, e1, e2, sigma, d, w] = deal(1000, 500, 0.80, 0.60, 1.712e-9, 1, 1);

% solve manually for some step size
[u, v]   = solveUV(0.5);
[I1, I2] = solveI(u, v);
[Q1, Q2] = solveQ(0.5);

% perform Richardson extrapolation
[Q, e]   = richardson(@solveQCombined, 1e-4)

function [u, v] = solveUV(h)

    global T1 T2 e1 e2 sigma d w

    % function handle
    F = @(x, y, d) 0.5 * 1 / (d^2 + (x - y)^2)^(3/2);

    % integral limits
    a = -0.5 * w;
    b =  0.5 * w;

    % step-size and division
    % h = 0.5;
    N = floor((b - a) / h);
    N = N + 1;

    % explanation:
    % the goal is to solve for u(x) and u(y)

    % in each equation, the integral is approximated using the Trapezoidal
    % method, resulting in a (N x 2N) matrix; for the u(x) equation,
    % each row corresponds to an equation for u(xi)
    
    % with u(x) = C1 + (1 - e1) * ∫ F(x, y, d) * v(y) dy as an example, for
    % each row, the difference equation is (approximately):
    % (1 - e1) * [2F(x1, y1) * v(y1) + F(x1, y2) * v(y2) ...] - u(x1) == - C1
    
    % two systems of linear equations, one per equation (u(x) and v(u)) is
    % setup and joined together as Az = b, where A = [Au; Av], b = [-C1; -C2]
    % and solved for z, which then contains [v[y1, y2 ..]; u[x1, x2, ..]]

    % Au = [      0.04  0.019046        -1          0]   [ v(y1) ]   [ -C1 ]
    %      [  0.019046      0.04         0         -1] * [ v(y2) ] = [ -C1 ]
    % Av = [        -1         0      0.08   0.038091]   [ u(x1) ]   [ -C2 ]
    %      [         0        -1  0.038091       0.08]   [ u(x2) ]   [ -C2 ]
    %      [-------------------------------------------|-----------|-------]
    %                             A                          z         b    

    % system of linear equations
    % Az = b

    Au = ones(N);
    Av = ones(N);

    for i = 1:N

        % u(xi) = C1 + (1 - e1) * ∫[-0.5w:h:0.5w] F(xi, yi, d) * v(yi) dyi
        % create row for Au[i, :] for a constant value xi and all yi
        yi = a + (i - 1) * h;
        Au(i, :) = (1 - e1) .* arrayfun(@(xi) F(xi, yi, d), a:h:b) .* trpz_weights(N, h);
        
        % v(yi) = C2 + (1 - e1) * ∫[-0.5w:h:0.5w] F(xi, yi, d) * u(xi) dxi
        % create row for Av[i, :] for a constant yi and all xi
        xi = a + (i - 1) * h;
        Av(i, :) = (1 - e2) .* arrayfun(@(yi) F(xi, yi, d), a:h:b) .* trpz_weights(N, h);

    end

    % combine A matrices diagonally
    A = blkdiag(Au, Av);

    % add -1 on the diagonals for -u(x) and -u(v) terms
    A = A + -1 * [ zeros(N) eye(N) ; eye(N) zeros(N)];

    % b: RHS vector
    C1 = (e1 * sigma * T1^4);
    C2 = (e2 * sigma * T2^4);
    b = [ones(N, 1) * -C1; ones(N, 1) * -C2];

    % solve for z using long division
    z = A \ b;

    % extract u, v vector
    v = z(1:N);
    u = z(N + 1: end);

end

function [I1, I2] = solveI(u, v)

    global T1 T2 e1 e2 sigma
    
    I1 = (u - e1 .* sigma .* T1.^4) ./ (1 - e1);
    I2 = (v - e2 .* sigma .* T2.^4) ./ (1 - e2);
end

function [Q1, Q2] = solveQ(h)

    [u, v]   = solveUV(h);
    [I1, I2] = solveI(u, v);
    
    Q1 = sum((u - I1) .* trpz_weights(length(u), h)');
    Q2 = sum((v - I2) .* trpz_weights(length(u), h)');

end

function Q = solveQCombined(h)

    [Q1, Q2] = solveQ(h);
    Q = [Q1; Q2];

end

function w = trpz_weights(N, h)

    % trapezoidal weights
    % h/2 * (f(x0) + 2f(x1) + ... 2(fx_N-1) + (fx_N))
    
    w = ones(1, N);
    w([2:end-1]) = 2;
    w = w .* (h / 2);
end

ode_bvp.m

% Second order Ordinary Differential Equations (ODE)
% Boundary Value Problem (BVP)

%%
%#ok<*NOPTS>
clear; close all; clc;

% config
global makeODE_prompt_gui
if ~exist('richardson')
    warning("To use richardson(), please include the .m file."); end

% example BVP problem
% y'' == 2x + sin(y') - cos(y)
% for x = [0; 2]
% y(0) = 0, y(2) = 1

% finite difference definitions
syms yprev yi ynext xi h
dy  = (ynext - yprev)/(2*h);
ddy = (yprev - 2*yi + ynext)/h^2;

% discrete ODE
eq = ddy == 2*xi + sin(dy) - cos(yi);
eq = eq - rhs(eq)

% check if ODE is on standard form
if ~has(eq, sym('0')) && (rhs(eq) ~= sym('0'))
    error("Specified ODE is NOT of standard form, e.g. y'' + cos(y') == 0"); end

% system boundaries
x  = [0 2]; % x boundaries: x = [a b]
y  = [0 1]; % boundary conditions: y = [y0 yN]
N  = 20;    % number of sub-divisions
Nx = N + 1; % number of ALL nodes (x-points)

if (exist('richardson') && questdlg("Solve using Richardson extrapolation?") == "Yes")

% use richardson extrapolation to solve y(1)
makeODE_prompt_gui = false;
sol = richardson(@(h) solveBVP(eq, 1, x, y, diff(x)/h), 1e-4);

else

% make system (returns function handles)
makeODE_prompt_gui = true;
[f, J] = makeODE(eq, x, y, N);

% column vector of initial guess using linear interpolation between 
% boundary conditions; spanning ALL nodes (not only internal nodes)
y0 = linspace(y(1), y(2), Nx)';

% solve and insert boundary conditiions
sol = newt(f, J, y0(2:end-1), 100);
sol = [ y(1); sol; y(2)];

% plot with interpolated x-axis
sol = [linspace(x(1), x(2), Nx)', sol];
plot(sol(:, 1), sol(:, 2))

% find value at y(1) (interpolate if necessary)
interp1(sol(:, 1), sol(:, 2), 1)

end

function [f, J] = makeODE(eq_sym, x, y, N)
    
    % constants
    a  = x(1);
    b  = x(2);
    y0 = y(1);
    yN = y(2);
    
    % define step size and redefine N
    % N is the number of total sub-divisions, redefined to the
    % number of internal nodes
    h = abs((b - a) / N);
    N = N - 1;

    % function vector
    F_sym = sym(zeros(N - 2, 1));
    
    % remove (== 0) if present
    if (has(eq_sym, sym('0')))
        eq_sym = lhs(eq_sym); end
       
    % first and last row
    syms yprev yi ynext
    F_sym(1) = subs(eq_sym, [yprev yi ynext], [y0 sym("y1") sym("y2")]);
    F_sym(N) = subs(eq_sym, [yprev yi ynext], [sym("y" + (N-1)) sym("y" + N) yN]);
    
    % intermediate rows
    for k = 2:N - 1
        F_sym(k) = subs(eq_sym, [yprev yi ynext], [sym("y" + (k-1)) sym("y" + k) sym("y" + (k+1))]);
    end
    
    % substitute h
    F_sym = subs(F_sym, sym('h'), h);
    
    % substitute xi
    for k = 1:N
        F_sym(k) = subs(F_sym(k), sym('xi'), a + k*h);
    end
    
    % Jacobian wrt. [y1, y2, .. yN]
    J_sym = jacobian(F_sym, sym('y', [1 N]));
    
    % function handles
    precision = 4;
    f = @(x) double(vpa(subs(F_sym, sym('y', [1 N]), x'), precision));
    J = @(x) double(vpa(subs(J_sym, sym('y', [1 N]), x'), precision));
    
    % show matrices
    global makeODE_prompt_gui
    if (makeODE_prompt_gui && questdlg("Show F and J matrix?") == "Yes")
        F_sym
        J_sym
    end
    
    % info
    disp("Created ODE system with step-size h = " + h + " with a total of " + N + " internal nodes.");

end

function [sol] = newt(f, J, x_init, max_iter)
    
    x_old = x_init;
    accuracy = 1e-6;
    
    disp("Solving system of non-linear equations using Newton-Raphson...");

    for i = 1:max_iter
        
        disp("i = " + i);
        
        dx = J(x_old) \ -f(x_old);
        x_new = x_old + dx;
        
        % update old estimate
        x_old = x_new;
        
        % check convergence
        if all(x_new ~= 0)
            ea = dx ./ x_new; end

        if max(abs(ea)) <= accuracy
            disp("Solution found!"), break; end
        
    end
    
    % check if solution was found
    if max(abs(ea)) > accuracy
            disp("No solution found!"); end
    
    % return solution
    sol = double(x_new);

end

function [sol] = solveBVP(eq_sym, desired_x, x, y, N)

    % make discrete ODE function handles
    [f, J] = makeODE(eq_sym, x, y, N);
    
    % initial guess for newt()
    y0 = linspace(y(1), y(2), N + 1)';
    
    % solve and insert boundary conditiions with interpolated x-axis
    sol = newt(f, J, y0(2:end-1), 100);
    sol = [ y(1); sol; y(2)];
    sol = [linspace(x(1), x(2), N + 1)', sol];
    
    % return solution for desired x value
    sol = interp1(sol(:, 1), sol(:, 2), desired_x);

end

ode_ivp.m

% First order Ordinary Differential Equations (ODE)
% Initial Value Problem (IVP)

%% analytical ODE
%#ok<*NOPTS>
clear; close all; clc;

% y’(x) = y(x) * cos(x + y(x)) with y(0)=1

syms y(x)

ode  = diff(y, x) == y(x) * cos(x + y(x))
cond = y(0) == 1

ysol(x) = dsolve(ode, cond)

%% system of ODE's (first order)
clc; close all; clear;

% system of first-order ODE's

% u'(x) = cos(-1 + x + u(x) + 3v(x))   | u(0) = 1
% v'(x) = -u(x)^2 + 2sin(v(x))         | v(0) = 0
 
% y'[]  = [ y0'(x) ] = [ cos(-1 + x + y0 + 3y1) ]
%         [ y1'(x) ] = [       -y0^2 + 2sin(y1) ]
 
% a[]   = [ 1 ]
%         [ 0 ]

% define ode, initial conditions, and range
ode   = @(x, y) [y(1) * cos(y(2)); -y(1)^3];
a     = [1, 0];
range = [0, 20];

% solve numerically
[x, y] = ode45(ode, range, a);

% print solution(s)
sol = y(end, :)

% plot solution(s)
plot(x, y, '-o')
legend("y[1]", "y[2]")

% jacobian wrt. y(x)
syms x y0 y1
J = jacobian([cos(-1 + x + y0 + 3 * y1); -y0^2 + 2 * sin(y1)], [y0, y1])

% eigen values
L = eig(J)

pde_parabolic_implicit.m

%#ok<*NOPTS>
clear; close all; clc;

% du/dt = a * d^2x/dx + f(x, t)
% f(x, t) = 10x * (1 - x) * exp(-t/10)
% a = 1

f = @(x, t) 10*x * (1 - x) * exp(-t/10);

% domain boundaries
x = [0 1];
y = [0 20]; % t
bounds = {@(x) x^4, @(y) 1, @(x) 0 , @(y) 0}; % TOP, RIGHT, BOTTOM, LEFT

% uniform step-size and number of sub-divisions
h = 0.2;
m = round(diff(y) / h) + 1;  % number of vertical nodes   (rows)
n = round(diff(x) / h) + 1 ; % number of horizontal nodes (cols)

% create matrix
Z = inf(m, n);

% apply boundary conditions, prioritize top/bottom
Z(:, 1)   = arrayfun(bounds{4}, y(1):h:y(2)); % LEFT
Z(:, end) = arrayfun(bounds{2}, y(1):h:y(2)); % RIGHT
% A(end, :) = arrayfun(bounds{3}, x(1):h:x(2)); % BOTTOM
Z(1, :)   = arrayfun(bounds{1}, x(1):h:x(2)); % TOP

% solve internal nodes, starting at initial row to solve next row
for row = 1:m - 1
	
	r = 1/(2*h);
	
	% -r * u(i-1, j+1) + (1 + 2r) * u(i, j+1) - r * u(i+1, j+1) = r * u(i-1, j) + (1 - 2*r) * u(ij) + r * u(i+1, j)

	% [ 1+2r    -r     0     0 ]   [ u22 ]   [ r * u11 + (1 - 2*r) * u21 + r * u31 - (-r * u12) ]
	% [   -r  1+2r    -r     0 ] * [ u32 ] = [ r * u21 + (1 - 2*r) * u31 + r * u41              ]
	% [    0    -r  1+2r    -r ]   [ u42 ]   [                        ...                       ]
	% [    0     0    -r  1+2r ]   [ u52 ]   [ r * u31 + (1 - 2*r) * u51 + r * u51 - (-r * u62) ]
	
	% compose A matrix
	A = [-r (1 + 2*r) -r];
	A = diag(ones(n - 3, 1) .* A(1), -1) + ...
	    diag(ones(n - 2, 1) .* A(2),  0) + ...
	    diag(ones(n - 3, 1) .* A(3),  1);
	
	% compose b matrix
	b = ones(n - 2, 1);
	for col = 2:n - 1
		
		xi = x(1) + (col - 1) * h;
		yi = y(1) + (row - 1) * h;
			   
		u_iprev_j  = Z(row, col - 1);
		u_ij       = Z(row, col);
		u_inext_j  = Z(row, col + 1);
		
		b(col - 1) = r*u_iprev_j + (1 - 2*r)*u_ij + r*u_inext_j + h*(f(xi, yi) + f(xi, yi + h));
		
	end
	
	% add boundary conditions
	b(1)   = b(1)   - (-r*Z(row, 1));
	b(end) = b(end) - (-r*Z(row, end));
	
	% solve
	z = A \ b;
	
	% insert solutions
	Z(row + 1, [2:end-1]) = z';
	
end

% interpolation
X = x(1):h:x(2);
Y = y(1):h:y(2);
u = @(x, y) interp2(X, Y, Z, x, y);

% surface plot
figure
surf(X, Y, Z, 'edgecolor', 'none')

% sliced plot
figure
imagesc(X, Y, Z)
colorbar

disp("u(0.5, 20) = " + u(0.5, 20))

pde_poisson.m

%#ok<*NOPTS>

function U = poisson_pde(f, u, lambda, x_range, y_range, bounds, h, xy, plot)

	% Poisson equation
	% -(uxx + uyy) + lambda * u(x, y) = f(x,y)

	% see example problem for usage
	% https://gist.github.com/martinandrovich/9eca71c2dd1d05bbb5ff9e9b4dd3317d#file-pde_poisson_demo-m

	lambda = 0;
	u = @(xi, yj) 0;
	f = @(xi, yj) (1 + xi + yj);

	% domain boundaries
	x0 = x_range(1);
	y0 = y_range(1);
	bound = containers.Map( ...
		{'TOP', 'RIGHT', 'BOTTOM', 'LEFT'}, [0 0 0 0]);

	% uniform step-size and number of sub-divisions
	n = floor(diff(x_range) / h); % number of horizontal subdivisions (cols)
	m = floor(diff(y_range) / h); % number of vertical subdivisions   (rows)

	disp("Solving Poisson PDE with step-size = " + h + "; using " + m + " x " + n + " sub-divisions.");

	% Z: solution matrix
	Z = zeros(m + 1, n + 1);

	% Z filled with boundary values
	Z(1, :) = bound('TOP');
	Z(:, end) = bound('RIGHT');
	Z(end, :) = bound('BOTTOM');
	Z(:, 1) = bound('LEFT');

	% system of linear equations
	% Aw = b

	% A: coefficient matrix [(Ny - 1)*(Nx - 1) x (Ny - 1)*(Nx - 1)]
	% A = zeros((n - 1)*(m - 1), (n - 1)*(m - 1));
	A = sparse((n - 1)*(m - 1), (n - 1)*(m - 1));

	% b: right hand size vector [(Ny - 1)*(Nx - 1) x 1]
	b = zeros((n - 1)*(m - 1), 1);

	% w vector (solutions)
	w = sym(b);

	% compose initial b and w vectors
	idx = 0;
	for i = 1:(m - 1)
		for j = 1:(n - 1)
			% increment index
			idx = idx + 1;

			% b vector
			xi = x0 + i * h;
			yj = y0 + j * h;
			b(idx) = h^2 * (f(xi, yj) - lambda * u(xi, yj));
		end
	end

	% compose A matrix and boundary conditions for b vector
	count = 0;
	for i = 1:(m - 1)*(n - 1)

		j = i;
		count = count + 1;

		% fill diagonal
		A(i, j) = 4 + h^2 * lambda;

		% u_(i, j-1) (left cell)
		if (count ~= 1)
			A(i, j - 1) = -1;
		else
			b(i) = b(i) + bound('LEFT');
		end

		% u_(i, j+1) (right cell)
		if (count ~= n - 1)
			A(i, j + 1) = -1;
		else
			b(i) = b(i) + bound('RIGHT');
		end

		% u_(i+1, j) (bottom cell)
		if(j + n - 1 <= (m - 1)*(n - 1))
			A(i, j + n - 1) = -1;
		else
			b(i) = b(i) + bound('BOTTOM');
		end

		% u_(i-1, j) (top cell)
		if (j - n + 1 > 0)
			A(i, j - n + 1) = -1;
		else
			b(i) = b(i) + bound('TOP');
		end

		% reached end of row
		if (count == n - 1)
			count = 0;
		end

	end

	% solve system
	w = A \ b;

	% insert into internal solutions into final stencil
	Z_internal = reshape(w, m - 1, n - 1);
	Z(2:m, 2:n) = Z_internal;

	% interpolation
	X = 0:h:diff(x_range);
	Y = 0:h:diff(y_range);

	if plot
		surf(X, Y , Z); end

	% interpolated solution at U(x, y)
	U = interp2(X, Y, Z, xy(1), xy(2));

end

pde_poisson_demo.m

clear; close all; clc;

% example PDE problem
% -(uxx + uyy) == 1 + x + y
%  (x, y) ∈  Ω where Ω = {(x, y) | 0 < x < 1, 0 < y < 1}
% u(x, y) ∈ dΩ where u(x, y) = 0

x_range = [0 1];
y_range = [0 1];
bounds  = [0 0 0 0]; % TOP, RIGHT, BOTTOM, LEFT
xy      = [0.5 0.5];
plot    = true;

lambda = 0;
u = @(xi, yj) 0;
f = @(xi, yj) (1 + xi + yj);

accuracy = 1e-5;
h0 = 0.8;

% solve for manual step-size and plot
U = pde_poisson(f, u, lambda, x_range, y_range, bounds, 0.1, xy, plot)

% solve using richardson extrapolation
[U0, e] = richardson(@(h) pde_poisson(f, u, lambda, x_range, y_range, bounds, h, xy, false), accuracy, h0)

richardson.m

%#ok<*NOPTS>
function [A0, e] = richardson(f, accuracy, h0, a)
    
    % function arguments validation
    % https://se.mathworks.com/help/matlab/matlab_prog/function-argument-validation-1.html
    
    arguments
        
        % function handle on the form @f(h), called with step-size
        % returning solution array sol[:]
        f (1,1) function_handle {mustBeValidFunctionHandle}
        
        % desired accuracy
        accuracy (1,1) double {mustBeNumeric, mustBeReal} = 1e-4
        
        % initial step size
        h0 (1,1) double {mustBePositive, mustBeReal} = 1.0
        
        % growth coefficient (alpha)
        a (1,1) double {mustBeInteger} = 2
        
        % function returns
        % A0[] = vector of optimal solutions
        % e[]  = vector of error for each solution
        % prints table when done
        
    end
    
    % internal varibales
    
    A = inf(3, 1);       % matrix of solution estimates; column n corresponds to y[n] (solution)
    h = [h0; 0; 0];      % vector of step-sizes, h0 < h1 < h2 && h1/h0 = h2/h1 = alpha
    e = inf(3, 1);       % vector of errors on A[h[1]]; error on solution n is e[n]
    k_est = inf(3, 1);   % vector of estimated orders k[]
    alpha_k = inf(3, 1); % vector of estimated alpha^k[]
    A0 = zeros(3, 1);    % vector of current most optimal solution estimates A0[]
    
    % output table
    T = array2table(zeros(0,7), "VariableNames", ["i", "h[0]", "A[h[0]]", "A0", "αᵏ", "~k", "e"]);
    
    % while while the absolute error is greater than required accuracy
    % or while k_est is not a finite value
    i = 0;
    num_warn = 0;
    while all(abs(e) > accuracy) || ~all(isfinite(k_est))
        
        try
            % compute solution(s)
            sol = f(h(1))';
        catch
            warning("Could not compute solution to ODE with step-size h = " + h(1) + "; increasing step-size.");
            num_warn = num_warn + 1;
            sol = NaN;
        end
        
        % too many warning; stop
        if num_warn > 10
            break; end
        
        % solution must be on column vector form
        if ~iscolumn(sol)
            sol = sol'; end
        
        % resize data if necessary
        % number of cols in 'A' must equal number of rows in 'sol'
        if size(A, 2) ~= size(sol, 1)
           A = inf(3, size(sol, 1));
           e = inf(size(sol));
           k_est = inf(size(sol));
           alpha_k = inf(size(sol));
           A0 = zeros(size(sol));
        end
        
        % propogate and insert new estimate(s)
		% the column-vector 'sol' is inserted as a row at A(1, :)
		% where A[1][n] is the most recent estimate of sol[n]
        A(3, :) = A(2, :);
        A(2, :) = A(1, :);
        A(1, :) = sol';
        
        % estimate coefficients
        alpha_k = ((A(3, :) - A(2, :)) ./  (A(2, :) - A(1, :)))';
        k_est = abs(log(alpha_k) / log(2));
        
        % estimated error
		e = 1 ./ (a.^k_est - 1)  .* (A(1, :) - A(2, :))';
        
        % richardson optimal estimate
		A0 = A(1, :)' + e;
        
        % insert table row      
        T = [T; {i, h(1), A(1, :), A0', alpha_k', k_est', e'}];
        
        % update step sizes
        i = i + 1;
        h(3) = h(2);
        h(2) = h(1);
        h(1) = h(1) / a;
        
        if h(1) <= 0
            error("Step-size is less than zero; aborting!"); end
             
    end
    
    % show table
    fprintf("\nRichardson extrapolation with desired accuracy of %.0d, an order of %d, and initial step-size of %0.1f:\n\n", accuracy, a, h0)
    disp(T)
    
end

function mustBeValidFunctionHandle(value)
    
    % checks for function_handle for richardson method
    % required: function [sol] = f(h)
    %           h = (1, 1), class(double)
    %           sol = (:, 1), class(double)
    
    if nargin(value) ~= 1
        error('Function handle must take exactly one argument as input.'); end
    
end