Please help me with syntax error on line 52

조회 수: 2 (최근 30일)
Saswat
Saswat 2024년 2월 29일
댓글: Walter Roberson 2024년 2월 29일
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)';

답변 (1개)

Walter Roberson
Walter Roberson 2024년 2월 29일
이동: Walter Roberson 2024년 2월 29일
% 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
Walter Roberson
Walter Roberson 2024년 2월 29일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Number Theory에 대해 자세히 알아보기

제품


릴리스

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by