Skip to content
Snippets Groups Projects
Commit 079781c7 authored by Nam Nguyen's avatar Nam Nguyen
Browse files

Add FD matlab code for RSHMHF

parents
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
README.md 0 → 100644
+ 10
− 0
View file @ 079781c7
# Code for solving harmonic map heat flow (HMHF) into 2-spheres
## Installation
[matlab](https://www.mathworks.com/products/matlab.html) are needed to run the code.
## Usage
- In `RSHMHF_FD`: Finite difference discretizations with [matlab](https://www.mathworks.com/products/matlab.html) for radially symmetric HMHF (RSHMHF) used in [arXiv:2503.03525](https://arxiv.org/abs/2503.03525)
RSHMHF_FD/RSHMHF_fd.m 0 → 100644
+ 271
− 0
View file @ 079781c7
close all; clear; clc;
%% Time parameters
bdf=1;
starttime = 0.0;
tstep = 1e-6;
endtime = 0.01;
n_tstep = 1e+4;
%% Grid for reference solution
N=999;
h_ref = 1/(N+1);
x_ref = linspace(0,1,N+2).';
%% choose initial condition and alpha
% case 1: u0 = pi * (1-x)*x; alpha = sqrt(2/pi)
% case 2: u0 = 9*pi * (1-x)*x; alpha = 1
choose_case = 2;
[u0_ref,alpha] = choose_example(choose_case,x_ref);
%% Numerical convergence test
do_time_conv = 0;
do_space_conv = 0;
%% Plot settings
plot_solution = 1;
do_energy = 1;
do_alpha = 1;
%% Compute reference solution
disp("Compute reference solution")
[sol_ref,energy,D_alpha_norm]= solveHMHF1D_FD(N,x_ref,u0_ref,bdf,...
starttime,tstep,endtime,n_tstep,alpha,do_energy,do_alpha);
disp("...done")
%% Plots
if plot_solution==1
figure
plot(x_ref,sol_ref,'LineWidth',5);
hold on;
plot(x_ref,u0_ref,'LineWidth',5);
hold off;
grid on;
xlabel('$x$','FontSize',40,'Interpreter','latex');
legend({'$\mathbf{u}^M$','$\mathbf{u}^0$'},"Interpreter","latex","FontSize",40);
xticks([0 0.2 0.4 0.6 0.8 1.0]);
yticks([0 2 4 6 8]);
ax = gca;
ax.FontSize=40;
ax.XAxis.TickLabels = {'$0$','$0.2$','$0.4$','$0.6$','$0.8$','$1.0$'};
ax.XAxis.TickLabelInterpreter='latex';
ax.YAxis.TickLabels = {'$0$','$2$','$4$','$6$','$8$'};
ax.YAxis.TickLabelInterpreter='latex';
end
if do_energy ==1
figure
time = starttime:tstep:endtime;
plot(time,energy,'LineWidth',5);
grid on;
legend({'$\mathcal{E}_h(\mathbf{u}^n)$'},"Interpreter","latex","FontSize",40);
xticks([0 0.002 0.004 0.006 0.008 0.01]);
yticks([90 100 110 120 130 140]);
xlabel('$t$','FontSize',40,'Interpreter','latex');
ax = gca;
ax.FontSize=40;
ax.XAxis.TickLabels = {'$0$','$0.002$','$0.004$','$0.006$','$0.008$','$0.01$'};
ax.XAxis.TickLabelInterpreter='latex';
ax.YAxis.TickLabels = {'$90$','$100$','$110$','$120$','$130$','$140$'};
ax.YAxis.TickLabelInterpreter='latex';
end
if do_alpha == 1
time = starttime:tstep:endtime;
figure
plot(time,D_alpha_norm,'LineWidth',5);
hold on;
fplot(@(x) D_alpha_norm(1),[starttime endtime],'LineWidth',5);
hold off;
grid on;
legend({'$\| \mathbf{u}^n_{\alpha} \|_{\infty}$','$\| \mathbf{u}^0_{\alpha} \|_{\infty}$'},"Interpreter","latex","FontSize",40);
xlabel('$t$','FontSize',40,'Interpreter','latex');
xticks([0 0.002 0.004 0.006 0.008 0.01]);
yticks([0 500 1000 1500 2000 2500]);
ax = gca;
ax.FontSize=40;
ax.XAxis.TickLabels = {'$0$','$0.002$','$0.004$','$0.006$','$0.008$','$0.01$'};
ax.XAxis.TickLabelInterpreter='latex';
ax.YAxis.TickLabels = {'$0$','$500$','$1000$','$1500$','$2000$','$2500$'};
ax.YAxis.TickLabelInterpreter='latex';
end
%% Convergence: time step
D = spdiags(x_ref,0,N+2,N+2);
if do_time_conv == 1
bdf=2;
tstep_arr = [1e-2, 1e-2/2, 1e-2/4, 1e-2/8, 1e-2/16];
err_arr = zeros(size(tstep_arr));
for i=1:length(tstep_arr)
tstep = tstep_arr(i);
disp("Solve");
[sol,~,~]= solveHMHF1D_FD(N,x_ref,u0_ref,bdf,starttime,tstep,...
endtime,n_tstep,alpha,do_energy,do_alpha);
disp("...done");
err_arr(i) = sqrt(h_ref*(D(2:N+1,2:N+1)*(sol_ref(2:N+1)-sol(2:N+1))).' *...
(sol_ref(2:N+1)-sol(2:N+1)));
end
eoc = zeros(length(err_arr)-1,1);
for i=1:length(err_arr)-1
eoc(i) = log(err_arr(i+1)/err_arr(i))/log(0.5);
end
err_arr
eoc
end
%% Convergence: mesh size
if do_space_conv==1
bdf=1;
N_arr = [2^2-1, 2^(3)-1, 2^(4)-1, 2^(5)-1, 2^(6)-1];
err_arr = zeros(size(N_arr));
for i=1:length(N_arr)
N = N_arr(i);
h = 1/(N+1);
x = linspace(0,1,N+2).';
[u0,~] = choose_example(choose_case,x);
D = spdiags(x,0,N+2,N+2);
D(1,1) = 1;
sol_ref_interp = interp1(x_ref,sol_ref,x,'cubic');
disp("Solve")
[sol,~,~]= solveHMHF1D_FD(N,x,u0,bdf,starttime,tstep,endtime,n_tstep,...
alpha,do_energy,do_alpha);
disp("...done")
err_arr(i) =sqrt( h*(D(2:N+1,2:N+1)*(sol_ref_interp(2:N+1)-sol(2:N+1))).' *...
(sol_ref_interp(2:N+1)-sol(2:N+1)));
end
eoc = zeros(length(err_arr)-1,1);
for i=1:length(err_arr)-1
eoc(i) = log(err_arr(i+1)/err_arr(i))/log(0.5);
end
err_arr
eoc
end
%% Solver
function [sol,energy,D_alpha_norm] = solveHMHF1D_FD(N,x,u0,bdf,starttime,tstep,endtime,ntstep,alpha,do_energy,do_alpha)
if alpha < 0
error("alpha needs to be positive!");
end
if bdf > 2
error("Only bdf 1 and 2 supported!");
elseif bdf < 0
error("Only bdf 1 and 2 supported!");
end
h = 1/(N+1);
D = spdiags(x,0,N+2,N+2);
D(1,1) = 1;
inv_D = inv(D);
inv_D_sq = inv_D.^2;
inv_D_alpha = inv_D.^alpha;
%%% second argument in diag: 1 = above diagonal; -1 = below diagonal
A =sparse(1/(h*h) *(diag(2*ones(1,N+2)) + diag(-ones(1,N+1),1) + diag(-ones(1,N+1),-1)));
B = sparse(1/(2*h) *(diag(zeros(1,N+2)) + diag(ones(1,N+1),1) + diag(-ones(1,N+1),-1)));
C=A-inv_D*B;
sol = u0;
sol_prev = u0;
if bdf==2
sol_prev_prev = u0;
end
t= starttime;
M = eye(2+N,2+N);
energy = [];
if do_energy ==1
energy = zeros(ntstep+1,1);
% initial energy
energy(1) = h*((D(2:N+1,2:N+1)*(C(2:N+1,2:N+1)*sol(2:N+1)) ).' * sol(2:N+1)+ ...
(inv_D(2:N+1,2:N+1)*F(sol(2:N+1))).' * F(sol(2:N+1)) );
end
D_alpha_norm = [];
if do_alpha ==1
D_alpha_norm = zeros(ntstep+1,1);
D_alpha_norm(1) = norm(inv_D_alpha(2:N+1,2:N+1) * sol(2:N+1),inf);
end
counter = 2;
if bdf==2
[sol,~,~] = solveHMHF1D_FD(N,x,u0,1,starttime,tstep,starttime+tstep,1,alpha,do_energy,do_alpha);
t = t+tstep;
sol_prev_prev(2:N+1) = sol_prev(2:N+1);
sol_prev(2:N+1) = sol(2:N+1);
if do_energy == 1
energy(counter) = h*((D(2:N+1,2:N+1)*(C(2:N+1,2:N+1)*sol(2:N+1)) ).' * sol(2:N+1)+ ...
(inv_D(2:N+1,2:N+1)*F(sol(2:N+1))).' * F(sol(2:N+1)) );
end
if do_alpha ==1
D_alpha_norm(counter) = norm(inv_D_alpha(2:N+1,2:N+1) * sol(2:N+1),inf);
end
counter = counter+1;
end
%%% initial condition needs to be zero at 0 for this to make sense:
%%% sin(2*0/(2*0) = 1 at x_1 = 0
g = sparse(eye(2+N,2+N));
while abs(t-endtime)>1e-12 && t<endtime
if bdf==1
g(2:N+2,2:N+2) = G(sol_prev(2:N+2));
M = eye(2+N,2+N)+tstep.*(C+g*inv_D_sq);
M(1,1) = 1;
M(1,2:N+2) = zeros(1,N+1);
M(N+2,N+2) = 1;
M(N+2,1:N+1) = zeros(1,N+1);
M = sparse(M);
elseif bdf==2
g(2:N+2,2:N+2) = G(2.*sol_prev(2:N+2)-sol_prev_prev(2:N+2));
M = 3.*eye(2+N,2+N)+2.*tstep.*(C+g*inv_D_sq);
M(1,1) = 3;
M(1,2:N+2) = zeros(1,N+1);
M(N+2,N+2) = 3;
M(N+2,1:N+1) = zeros(1,N+1);
M = sparse(M);
end
if bdf==1
sol = M\sol_prev;
elseif bdf==2
sol = M\(4.*sol_prev-sol_prev_prev);
sol_prev_prev(2:N+1) = sol_prev(2:N+1);
end
sol_prev(2:N+1) = sol(2:N+1);
if do_energy == 1
energy(counter) = h*((D(2:N+1,2:N+1)*(C(2:N+1,2:N+1)*sol(2:N+1)) ).' * sol(2:N+1)+ ...
(inv_D(2:N+1,2:N+1)*F(sol(2:N+1))).' * F(sol(2:N+1)) );
end
if do_alpha==1
D_alpha_norm(counter) = norm(inv_D_alpha(2:N+1,2:N+1) * sol(2:N+1),inf);
end
counter = counter+1;
t=t+tstep;
%disp(t)
end
end
function [y] = G(u)
y = sparse(diag(sin(2.*u)./(2.*u)));
end
function [y] = F(u)
y = sin(u);
end
function [u0,alpha] = choose_example(which_case,x)
if which_case ==1
u0 = pi * (1-x).*x;
alpha = sqrt(2/pi);
elseif which_case==2
u0 = 9*pi * (1-x).*x;
%alpha = sqrt(2/(9*pi));
alpha = 1;
else
error('No case for initial condition defined!');
end
end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment