Effacer les filtres
Effacer les filtres

problem onn numerical solution

1 vue (au cours des 30 derniers jours)
ABDO
ABDO le 26 Avr 2024
Commenté : ABDO le 26 Avr 2024
= 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
  5 commentaires
Torsten
Torsten le 26 Avr 2024
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
ABDO le 26 Avr 2024
Smoothing length h is a length that allows the support of a cut-off core to be fixed W

Connectez-vous pour commenter.

Réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by