You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
I wrote code on MATLAB to solve Stokes equations using Bramble and Pasciak method but I got strange values of error for pressure p
3 views (last 30 days)
Show older comments
%Define problem parameters
mu = 1.0; % viscosity
%f = [1; 0; 1; 2]; % forcing term for velocity
%g = [2; 1]; % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
Unrecognized function or variable 'squaremesh'.
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B'; % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1); % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1)) 2*(eqn.A)\eqn.B_transpose; eqn.B 2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A zeros(size(transpose(eqn.B))); zeros(size(eqn.B)) eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
if H1ErrUCg(i) == 0 && i == 1
H1ErrUCg(i) = max(H1ErrUCg);
elseif H1ErrUCg(i) == 0
H1ErrUCg(i) = H1ErrUCg(i-1);
end
end
for i = 1:iterCg
if L2ErrPCg(i) == 0 && i == 1
L2ErrPCg(i) = max(L2ErrPCg);
elseif L2ErrPCg(i) == 0
L2ErrPCg(i) = L2ErrPCg(i-1);
end
end
for i = 1:iterCg
if L2ErrGradUandPCg(i) == 0 && i == 1
L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
elseif L2ErrGradUandPCg(i) == 0
L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) * C * x;
end
17 Comments
Mostafa Kassoum
on 12 Jan 2024
Edited: Walter Roberson
on 12 Jan 2024
The sequare mesh comes from the following link;
Mostafa Kassoum
on 12 Jan 2024
Edited: Walter Roberson
on 12 Jan 2024
Also to test the code you need the following code:
function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
z_old = x;
p_old = F - M*z_old;
r_old=p_old;
for i = 1:maxiter
alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
z_new = z_old + alpha_i*p_old;
r_new = F - M*z_new;
norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
L2ErrP(i) = norm(p - z_new(length(u)+1:end));
L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
%L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
relres(i) = norm_r;
fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
if (norm_r)< tol
break
end
beta_i = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
p_new = r_new - beta_i*p_old;
% Update
z_old = z_new;
r_old = r_new;
p_old = p_new;
end
end
Walter Roberson
on 12 Jan 2024
The given archive does not include ExtractBoundaryEdges or StokesModelMustafaNew
%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";
tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");
mkdir(tn);
addpath(tn)
websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
%done bringing in external tools
%Define problem parameters
mu = 1.0; % viscosity
%f = [1; 0; 1; 2]; % forcing term for velocity
%g = [2; 1]; % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1588886/image.png)
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
Unrecognized function or variable 'ExtractBoundaryEdges'.
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B'; % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1); % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1)) 2*(eqn.A)\eqn.B_transpose; eqn.B 2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A zeros(size(transpose(eqn.B))); zeros(size(eqn.B)) eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
if H1ErrUCg(i) == 0 && i == 1
H1ErrUCg(i) = max(H1ErrUCg);
elseif H1ErrUCg(i) == 0
H1ErrUCg(i) = H1ErrUCg(i-1);
end
end
for i = 1:iterCg
if L2ErrPCg(i) == 0 && i == 1
L2ErrPCg(i) = max(L2ErrPCg);
elseif L2ErrPCg(i) == 0
L2ErrPCg(i) = L2ErrPCg(i-1);
end
end
for i = 1:iterCg
if L2ErrGradUandPCg(i) == 0 && i == 1
L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
elseif L2ErrGradUandPCg(i) == 0
L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) * C * x;
end
function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
z_old = x;
p_old = F - M*z_old;
r_old=p_old;
for i = 1:maxiter
alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
z_new = z_old + alpha_i*p_old;
r_new = F - M*z_new;
norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
L2ErrP(i) = norm(p - z_new(length(u)+1:end));
L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
%L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
relres(i) = norm_r;
fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
if (norm_r)< tol
break
end
beta_i = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
p_new = r_new - beta_i*p_old;
% Update
z_old = z_new;
r_old = r_new;
p_old = p_new;
end
end
Mostafa Kassoum
on 12 Jan 2024
Edited: Walter Roberson
on 12 Jan 2024
The following code is the StokesModelMustafaNew
function pde = StokesModelMustafaNew
%% Custom Stokes equation model
%
pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
% exact velocity
function z = exactu(p)
x = p(:,1);
y = p(:,2);
z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
end
% exact pressure
function z = exactp(p)
x = p(:,1); y = p(:,2);
z = 2/3-x.^2-y.^2;
end
% Dirichlet boundary condition of velocity
function z = g_D(p)
z = exactu(p);
end
function z = g_N(p)
z(:,1) = 2*ones(size(p,1),1);
z(:,2) = zeros(size(p,1),1);
end
end
Mostafa Kassoum
on 12 Jan 2024
I think now you can run the code and see why is the L^2 error of pressure p big value.
Walter Roberson
on 12 Jan 2024
Still need ExtractBoundaryEdges
%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";
tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");
mkdir(tn);
addpath(tn)
websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
%done bringing in external tools
%Define problem parameters
mu = 1.0; % viscosity
%f = [1; 0; 1; 2]; % forcing term for velocity
%g = [2; 1]; % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1589536/image.png)
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
Unrecognized function or variable 'ExtractBoundaryEdges'.
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B'; % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1); % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1)) 2*(eqn.A)\eqn.B_transpose; eqn.B 2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A zeros(size(transpose(eqn.B))); zeros(size(eqn.B)) eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
if H1ErrUCg(i) == 0 && i == 1
H1ErrUCg(i) = max(H1ErrUCg);
elseif H1ErrUCg(i) == 0
H1ErrUCg(i) = H1ErrUCg(i-1);
end
end
for i = 1:iterCg
if L2ErrPCg(i) == 0 && i == 1
L2ErrPCg(i) = max(L2ErrPCg);
elseif L2ErrPCg(i) == 0
L2ErrPCg(i) = L2ErrPCg(i-1);
end
end
for i = 1:iterCg
if L2ErrGradUandPCg(i) == 0 && i == 1
L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
elseif L2ErrGradUandPCg(i) == 0
L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) * C * x;
end
function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
z_old = x;
p_old = F - M*z_old;
r_old=p_old;
for i = 1:maxiter
alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
z_new = z_old + alpha_i*p_old;
r_new = F - M*z_new;
norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
L2ErrP(i) = norm(p - z_new(length(u)+1:end));
L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
%L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
relres(i) = norm_r;
fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
if (norm_r)< tol
break
end
beta_i = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
p_new = r_new - beta_i*p_old;
% Update
z_old = z_new;
r_old = r_new;
p_old = p_new;
end
end
function pde = StokesModelMustafaNew
%% Custom Stokes equation model
%
pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
% exact velocity
function z = exactu(p)
x = p(:,1);
y = p(:,2);
z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
end
% exact pressure
function z = exactp(p)
x = p(:,1); y = p(:,2);
z = 2/3-x.^2-y.^2;
end
% Dirichlet boundary condition of velocity
function z = g_D(p)
z = exactu(p);
end
function z = g_N(p)
z(:,1) = 2*ones(size(p,1),1);
z(:,2) = zeros(size(p,1),1);
end
end
Mostafa Kassoum
on 12 Jan 2024
The following code is ExtractBoundaryEdges.
function [Dirichlet,bdFlag] = ExtractBoundaryEdges(elem, node)
bdFlag = setboundary(node,elem,'Dirichlet','(x==-1) | (x==1) | (y==1) | (y==-1)');%,'Neumann', '(y==1)');
% bdFlag = setboundary(node,elem,'Neumann','(y==1)');
% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
Dirichlet = allEdge((bdFlag(:) == 1),:);
end
Walter Roberson
on 12 Jan 2024
Need StokesisoP2P0
%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";
setb_url = "https://raw.githubusercontent.com/lyc102/ifem/master/fem/setboundary.m";
tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");
setb_localname = fullfile(tn, "setboundary.m");
mkdir(tn);
addpath(tn)
websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
websave(setb_localname, setb_url);
%done bringing in external tools
%Define problem parameters
mu = 1.0; % viscosity
%f = [1; 0; 1; 2]; % forcing term for velocity
%g = [2; 1]; % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1589696/image.png)
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
Unrecognized function or variable 'StokesisoP2P0'.
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B'; % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1); % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1)) 2*(eqn.A)\eqn.B_transpose; eqn.B 2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A zeros(size(transpose(eqn.B))); zeros(size(eqn.B)) eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
if H1ErrUCg(i) == 0 && i == 1
H1ErrUCg(i) = max(H1ErrUCg);
elseif H1ErrUCg(i) == 0
H1ErrUCg(i) = H1ErrUCg(i-1);
end
end
for i = 1:iterCg
if L2ErrPCg(i) == 0 && i == 1
L2ErrPCg(i) = max(L2ErrPCg);
elseif L2ErrPCg(i) == 0
L2ErrPCg(i) = L2ErrPCg(i-1);
end
end
for i = 1:iterCg
if L2ErrGradUandPCg(i) == 0 && i == 1
L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
elseif L2ErrGradUandPCg(i) == 0
L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) * C * x;
end
function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
z_old = x;
p_old = F - M*z_old;
r_old=p_old;
for i = 1:maxiter
alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
z_new = z_old + alpha_i*p_old;
r_new = F - M*z_new;
norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
L2ErrP(i) = norm(p - z_new(length(u)+1:end));
L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
%L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
relres(i) = norm_r;
fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
if (norm_r)< tol
break
end
beta_i = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
p_new = r_new - beta_i*p_old;
% Update
z_old = z_new;
r_old = r_new;
p_old = p_new;
end
end
function pde = StokesModelMustafaNew
%% Custom Stokes equation model
%
pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
% exact velocity
function z = exactu(p)
x = p(:,1);
y = p(:,2);
z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
end
% exact pressure
function z = exactp(p)
x = p(:,1); y = p(:,2);
z = 2/3-x.^2-y.^2;
end
% Dirichlet boundary condition of velocity
function z = g_D(p)
z = exactu(p);
end
function z = g_N(p)
z(:,1) = 2*ones(size(p,1),1);
z(:,2) = zeros(size(p,1),1);
end
end
function [Dirichlet,bdFlag] = ExtractBoundaryEdges(elem, node)
bdFlag = setboundary(node,elem,'Dirichlet','(x==-1) | (x==1) | (y==1) | (y==-1)');%,'Neumann', '(y==1)');
% bdFlag = setboundary(node,elem,'Neumann','(y==1)');
% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
Dirichlet = allEdge((bdFlag(:) == 1),:);
end
Mostafa Kassoum
on 13 Jan 2024
You need download the codes in the following link on MATLAB to get what you need.
Walter Roberson
on 13 Jan 2024
%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";
setb_url = "https://raw.githubusercontent.com/lyc102/ifem/master/fem/setboundary.m";
S2P0_url = "https://raw.githubusercontent.com/lyc102/ifem/master/equation/StokesisoP2P0.m";
dof_url = "https://raw.githubusercontent.com/lyc102/ifem/master/dof/dof3P2.m";
tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");
setb_localname = fullfile(tn, "setboundary.m");
S2P0_localname = fullfile(tn, "StokesisoP2P0.m");
dof3_localname = fullfile(tn, "dof3P2.m");
dof_localname = fullfile(tn, "dofP2.m");
mkdir(tn);
addpath(tn)
websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
websave(setb_localname, setb_url);
websave(S2P0_localname, S2P0_url);
websave(dof_localname, dof_url);
websave(dof3_localname, dof_url);
%done bringing in external tools
%Define problem parameters
mu = 1.0; % viscosity
%f = [1; 0; 1; 2]; % forcing term for velocity
%g = [2; 1]; % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1589746/image.png)
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000;
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
Index in position 2 exceeds array bounds. Index must not exceed 3.
Error in dofP2 (line 17)
totalEdge = uint32([elem(:,[1 2]); elem(:,[1 3]); elem(:,[1 4]); ...
Error in StokesisoP2P0 (line 25)
[tempvar,edgeC] = dofP2(elemC);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B'; % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1); % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1)) 2*(eqn.A)\eqn.B_transpose; eqn.B 2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A zeros(size(transpose(eqn.B))); zeros(size(eqn.B)) eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
if H1ErrUCg(i) == 0 && i == 1
H1ErrUCg(i) = max(H1ErrUCg);
elseif H1ErrUCg(i) == 0
H1ErrUCg(i) = H1ErrUCg(i-1);
end
end
for i = 1:iterCg
if L2ErrPCg(i) == 0 && i == 1
L2ErrPCg(i) = max(L2ErrPCg);
elseif L2ErrPCg(i) == 0
L2ErrPCg(i) = L2ErrPCg(i-1);
end
end
for i = 1:iterCg
if L2ErrGradUandPCg(i) == 0 && i == 1
L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
elseif L2ErrGradUandPCg(i) == 0
L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) * C * x;
end
function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
z_old = x;
p_old = F - M*z_old;
r_old=p_old;
for i = 1:maxiter
alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
z_new = z_old + alpha_i*p_old;
r_new = F - M*z_new;
norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
L2ErrP(i) = norm(p - z_new(length(u)+1:end));
L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
%L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
relres(i) = norm_r;
fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
if (norm_r)< tol
break
end
beta_i = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
p_new = r_new - beta_i*p_old;
% Update
z_old = z_new;
r_old = r_new;
p_old = p_new;
end
end
function pde = StokesModelMustafaNew
%% Custom Stokes equation model
%
pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
% exact velocity
function z = exactu(p)
x = p(:,1);
y = p(:,2);
z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
end
% exact pressure
function z = exactp(p)
x = p(:,1); y = p(:,2);
z = 2/3-x.^2-y.^2;
end
% Dirichlet boundary condition of velocity
function z = g_D(p)
z = exactu(p);
end
function z = g_N(p)
z(:,1) = 2*ones(size(p,1),1);
z(:,2) = zeros(size(p,1),1);
end
end
function [Dirichlet,bdFlag] = ExtractBoundaryEdges(elem, node)
bdFlag = setboundary(node,elem,'Dirichlet','(x==-1) | (x==1) | (y==1) | (y==-1)');%,'Neumann', '(y==1)');
% bdFlag = setboundary(node,elem,'Neumann','(y==1)');
% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
Dirichlet = allEdge((bdFlag(:) == 1),:);
end
Mostafa Kassoum
on 13 Jan 2024
It worked for me anf I got the results for everything and they are nice results except the error of pressure P, they are so big. I think there is a mistake but I didn't got it.
Mostafa Kassoum
on 13 Jan 2024
Foe example I got the following results:
H1 error for velocity u are: 0.129922159092270 0.138027304457401
L2 error for pressure p are: 12.6454430756918 12.6450517755845
I need know why are L2 error p big
Walter Roberson
on 13 Jan 2024
You can copy my entire post to your system and test it locally.
I suspect that something you have installed does not match the latest code at https://github.com/lyc102/ifem
Mostafa Kassoum
on 23 Jan 2024
I copied your entire post and runned it but I got also a big value of L2 error 11.3255614294248 11.3251245255107
As you can see they are almost equal to the values I obtained previously. Where do you think the mistake is?
Answers (0)
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)