Skip to content

Instantly share code, notes, and snippets.

@martinandrovich

martinandrovich/mandatory_2.m

Last active June 8, 2020 19:16
Show Gist options
  • Save martinandrovich/9eca71c2dd1d05bbb5ff9e9b4dd3317d to your computer and use it in GitHub Desktop.
Save martinandrovich/9eca71c2dd1d05bbb5ff9e9b4dd3317d to your computer and use it in GitHub Desktop.
Download ZIP
ODE and PDE implementations in MATLAB
Raw
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
Raw
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
Raw
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)
Raw
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))
Raw
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
Raw
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)
Raw
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment