Borrar filtros
Borrar filtros

Please help me with syntax error on line 52

8 visualizaciones (últimos 30 días)
Saswat
Saswat el 29 de Feb. de 2024
Comentada: Walter Roberson el 29 de Feb. de 2024
function [e1, e2] = spand_hw7(mesh, d, f, dfdx)
% Check if input arguments are provided
if nargin < 4
error('Insufficient input arguments. Please provide all required inputs.');
end
% Initialize error norms
e1 = 0;
e2 = 0;
% Define Gauss quadrature points and weights
x_quad = [-sqrt(5 + 2*sqrt(10/7))/3, 0, sqrt(5 + 2*sqrt(10/7))/3];
w_quad = [5/9, 8/9, 5/9];
% Loop over elements
for el = 1:size(mesh.conn, 2)
% Element nodes
nodes = mesh.conn(:, el);
% Element coordinates
xe = mesh.x(nodes);
% Element displacements
de = d(nodes);
% Jacobian matrix
J = (xe(end) - xe(1)) / 2;
% Initialize element error norms
e1_el = 0;
e2_el = 0;
% Loop over quadrature points
for q = 1:length(x_quad)
% Quadrature point
xq = (xe(end) + xe(1)) / 2 + x_quad(q) * J;
% Interpolated solution at quadrature point
u_interp = sum(N(xq, xe) .* de');
% Exact solution and its derivative at quadrature point
u_exact = f(xq);
du_exact = dfdx(xq);
% Compute L2 error norm
e1_el = e1_el + w_quad(q) * (u_exact - u_interp)^2;
% Compute energy error norm
e2_el = e2_el + w_quad(q) * sum((du_exact - de' .* dNdx(xq, xe)).^2);
end
% Add element error norms to total error norms
e1 = e1 + e1_el * J;
e2 = e2 + e2_el * J;
end
% Normalize error norms
e1 = sqrt(e1);
e2 = sqrt(e2);
end
% Shape functions
function N_val = N(x, xe)
N_val = [(x - xe(2)) .* (x - xe(3)) / ((xe(1) - xe(2)) * (xe(1) - xe(3)));
(x - xe(1)) .* (x - xe(3)) / ((xe(2) - xe(1)) * (xe(2) - xe(3)));
(x - xe(1)) .* (x - xe(2)) / ((xe(3) - xe(1)) * (xe(3) - xe(2)))];
end
% Derivative of shape functions
function dNdx_val = dNdx(x, xe)
dNdx_val = [(2*x - xe(2) - xe(3)) / ((xe(1) - xe(2)) * (xe(1) - xe(3)));
(2*x - xe(1) - xe(3)) / ((xe(2) - xe(1)) * (xe(2) - xe(3)));
(2*x - xe(1) - xe(2)) / ((xe(3) - xe(1)) * (xe(3) - xe(2)))];
end
Syntax error at line 52: Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To operate on each element of the matrix individually, use POWER (.^) for elementwise power.
% Example inputs before calling the function
ne = 10;
nn = 2*ne + 1;
mesh.x = linspace(0, 6, nn);
mesh.conn = [1:2:nn-2; 2:2:nn-1; 3:2:nn];
f = @(x) 1.0./((x-3).^2+0.2);
dfdx = @(x) -2.0*(x-3.0) ./ ((x-3).^2 + 0.2).^2;
d = f(mesh.x)';

Respuestas (1)

Walter Roberson
Walter Roberson el 29 de Feb. de 2024
Movida: Walter Roberson el 29 de Feb. de 2024
% Example inputs before calling the function
ne = 10;
nn = 2*ne + 1;
mesh.x = linspace(0, 6, nn);
mesh.conn = [1:2:nn-2; 2:2:nn-1; 3:2:nn];
f = @(x) 1.0./((x-3).^2+0.2);
dfdx = @(x) -2.0*(x-3.0) ./ ((x-3).^2 + 0.2).^2;
d = f(mesh.x)';
[E1,E2] = spand_hw7(mesh, d, f, dfdx)
de = 3×1
0.1087 0.1335 0.1678
Name Size Bytes Class Attributes Nxq 3x1 24 double Name Size Bytes Class Attributes u_exact 1x1 8 double u_interp 1x3 24 double
Error using ^
Incorrect dimensions for raising a matrix to a power. Check that the matrix is square and the power is a scalar. To operate on each element of the matrix individually, use POWER (.^) for elementwise power.

Error in solution>spand_hw7 (line 61)
e1_el = e1_el + w_quad(q) * (u_exact - u_interp)^2;
function [e1, e2] = spand_hw7(mesh, d, f, dfdx)
% Check if input arguments are provided
if nargin < 4
error('Insufficient input arguments. Please provide all required inputs.');
end
% Initialize error norms
e1 = 0;
e2 = 0;
% Define Gauss quadrature points and weights
x_quad = [-sqrt(5 + 2*sqrt(10/7))/3, 0, sqrt(5 + 2*sqrt(10/7))/3];
w_quad = [5/9, 8/9, 5/9];
% Loop over elements
for el = 1:size(mesh.conn, 2)
% Element nodes
nodes = mesh.conn(:, el);
% Element coordinates
xe = mesh.x(nodes);
% Element displacements
de = d(nodes)
% Jacobian matrix
J = (xe(end) - xe(1)) / 2;
% Initialize element error norms
e1_el = 0;
e2_el = 0;
% Loop over quadrature points
for q = 1:length(x_quad)
% Quadrature point
xq = (xe(end) + xe(1)) / 2 + x_quad(q) * J;
% Interpolated solution at quadrature point
Nxq = N(xq, xe);
whos Nxq
u_interp = sum(N(xq, xe) .* de');
% Exact solution and its derivative at quadrature point
u_exact = f(xq);
du_exact = dfdx(xq);
whos u_exact u_interp
% Compute L2 error norm
e1_el = e1_el + w_quad(q) * (u_exact - u_interp)^2;
% Compute energy error norm
e2_el = e2_el + w_quad(q) * sum((du_exact - de' .* dNdx(xq, xe)).^2);
end
% Add element error norms to total error norms
e1 = e1 + e1_el * J;
e2 = e2 + e2_el * J;
end
% Normalize error norms
e1 = sqrt(e1);
e2 = sqrt(e2);
end
% Shape functions
function N_val = N(x, xe)
N_val = [(x - xe(2)) .* (x - xe(3)) / ((xe(1) - xe(2)) * (xe(1) - xe(3)));
(x - xe(1)) .* (x - xe(3)) / ((xe(2) - xe(1)) * (xe(2) - xe(3)));
(x - xe(1)) .* (x - xe(2)) / ((xe(3) - xe(1)) * (xe(3) - xe(2)))];
end
% Derivative of shape functions
function dNdx_val = dNdx(x, xe)
dNdx_val = [(2*x - xe(2) - xe(3)) / ((xe(1) - xe(2)) * (xe(1) - xe(3)));
(2*x - xe(1) - xe(3)) / ((xe(2) - xe(1)) * (xe(2) - xe(3)));
(2*x - xe(1) - xe(2)) / ((xe(3) - xe(1)) * (xe(3) - xe(2)))];
end
  1 comentario
Walter Roberson
Walter Roberson el 29 de Feb. de 2024
u_interp = sum(N(xq, xe) .* de');
N(xq, xe) is 3 x 1
de is 3 x 1
de' is 1 x 3
3 x 1 .* 1 x 3 gives 3 x 3
sum(3 x 3) gives 1 x 3
So u_interp comes out 1 x 3.

Iniciar sesión para comentar.

Productos


Versión

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by