problem onn numerical solution

= 0.3; % Rayon de lissage
dx=0.3;
N =floor(1/dx);
x =0:dx:(pi/sqrt(3)); % Positions des particules
dt = 1/4;
t = 0:dt:1;
M = length(t) - 1;
% Initialisation de la solution approchée et de la matrice A
A = zeros(N, N);
b=zeros(N,1);
uapp=zeros(N,1);
%sol et cond exacte
u_exact = @(x) cos(sqrt(3)*x);
% Construction de la matrice A et du vecteur b
for i = 1:N
    for j =N
        r_ij = abs(x(i) - x(j));
        if r_ij <= h         
            d2W_dx2 = second_derivative_monaghan_kernel(r_ij, h);
            W_value = monaghan_kernel(r_ij, h);
            A(i, j) = (u_exact(x(j))- u_exact(x(i))) * d2W_dx2;
            b(i) = -3 * u_exact(x(j)) * W_value ;
        end
    end
end 
u_app(1)=1;
u_app(N)=u_exact(pi/sqrt(3));
uapp = b\A; 
uapp = [0; uapp];
% Affichage des résultats
plot(x , uapp, 'r');
hold on;
plot(x, u_exact(x), 'b--'); 
xlabel('x');
ylabel('u');
legend('Solution approchée (SPH)', 'Solution exacte');
title('Solution approchée et exacte par SPH');
grid on;
% Fonctions du noyau de Monaghan
function W_value =monaghan_kernel(r_ij, h)
C= (1/ h);
if r_ij < h
    W_value =C*((2 / 3) - (r_ij / h)^2 +(1/2 *(r_ij / h)^3));
elseif r_ij<2*h
    W_value =C*((1/6)*(2-(r_ij / h))^3);      
else
    W_value=0;
end 
end
function d2W_dx2 = second_derivative_monaghan_kernel(r_ij, h)
C = (1 / h); % Corrected the denominator
if r_ij < h
    d2W_dx2 = -2 * C * (1 - (3 * (r_ij / h)));
elseif r_ij < 2*h
    d2W_dx2 = 2 * C * (2 - (r_ij / h)); % Removed extra characters
else
    d2W_dx2 = 0;
end
end

I have a problem simulataing a numerical solution for this program

Ineed some help urgent

[Steven Lord]

What problems are you experiencing?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

[Walter Roberson]

When you posted the code, you accidentally left out the left hand side of the first assignment.

[ABDO]

the numerical solution must be closer than the exact solution but I'm having a problem displaying the solution

maybe one of the terms is inacalculable

[Torsten]

It would help if you explained the problem you are trying to solve.

Further - as @Walter Roberson mentionned : the first line of your code is incomplete. Did you mean

h = 0.3; % Rayon de lissage

?

[ABDO]

Smoothing length h is a length that allows the support of a cut-off core to be fixed W

Answer 1

[Vinay]

Hii ABDO,

The Smoothed Particle Hydrodynamics (SPH) method calculates the approximate values by solving the linear equation. The error arises in the code due to the following issues.

• Loop Mistake: In the nested loop, the variable iterates for ‘j’ = ‘N’, which should be for ‘j’ = 1:N. This means the loop is only iterating over the last element of ‘x’, causing an error.

% Construction of matrix A and vector b
for i = 1:N
    for j = 1:N % Corrected loop range
        r_ij = abs(x(i) - x(j));
        if r_ij <= h
            d2W_dx2 = second_derivative_monaghan_kernel(r_ij, h);
            W_value = monaghan_kernel(r_ij, h);
            % Update A(i, j) based on the interaction
            A(i, j) = d2W_dx2;
            % Accumulate b(i) with contributions from j
            b(i) = b(i) - 3 * u_exact(x(j)) * W_value;
        end
    end
end 

• Matrix Solver: The calculation ‘uapp’ = b\A is incorrect because ‘b’ is a column vector and ‘A’ is a square matrix. The solver solves the ‘A \ b’ to find ‘uapp’.

% Solve system
uapp = A \ b; % Solve the linear system

• Incorrect dimension: The dimension of the ‘x’ and the ‘uapp’ are not compatible, the dimension of ‘uapp’ should be same as that of ‘x’.

x = 0:dx:(pi/sqrt(3)); % Positions of particles
N = length(x); % Update N to match the length of x