C'est une partie d'un problème PDE que j'essaie de résoudre numériquement. Je sais que mes calculs pour les entrées sont corrects, mais j'ai du mal à assembler correctement la matrice dans Matlab. Les premières lignes sont correctes, mais les deux dernières lignes ne sont pas correctes car il me manque les éléments a_{11,10}, a_{10,10} et a_{10,11}. Ils me manquent car ma boucle n'y accède jamais. Par exemple a_{10,10} nécessite i = 10, qui est un nombre pair, cependant mon M = 11 dans ce cas et je ne passe que de 1 à 9. Si je fais passer ma boucle for de 1:11, je sortir de l'erreur des limites puisque j'utilise les indecies i+1 et i+2.

Il en va de même pour le vecteur b. Comment puis-je remplir tous mes b_i si je ne passe qu'à 9 mais que le vecteur fait 11 éléments ?

Une suggestion sur la façon de résoudre ce problème ? J'ai essayé de créer des boucles for séparées pour les cas pair et impair, mais j'ai toujours le même problème d'entrées manquantes en raison de l'utilisation des indecies i+1 et i+2. Ou dois-je saisir manuellement (code en dur) ces indices en dehors de la boucle foor ?

Voici ma boucle:


h = 1/10; 
x = 0:h:1;
M = length(x);  

% Assemble stiffness matrix A and load vector b.

A = zeros(M,M);
b = zeros(M,1);

for i = 1:M-2 
    if rem(i,2) ~= 0                   % if i is odd
        A(i,i) = A(i,i) + 1/(3*h);

        b(i) = b(i) + f(x(i)) * (-2)*h/3;

    elseif rem(i,2) == 0               % if i is even
        A(i,i) = A(i,i) + 3/(3*h);
        A(i+2,i) = A(i+2,i) - 8/(3*h);
        A(i,i+2) = A(i,i+2) - 8/(3*h);

        b(i) = b(i) + f(x(i)) * 8*h/3;
    end

    A(i+1,i) = A(i+1,i) + 2/(3*h);
    A(i,i+1) = A(i,i+1) - 3/(3*h);
end

% Enforce the boundary conditions

A(1,1) = 1.e4;
A(end,end) = 1.e4;

0
Parseval 9 févr. 2020 à 22:54

1 réponse

Meilleure réponse

Je pense que vous pouvez simplement remplir la matrice A avec deux lignes et colonnes supplémentaires, et les supprimer plus tard :

Regardez la version modifiée suivante de votre code :

h = 1/10; 
x = 0:h:1;
M = length(x);  

% Assemble stiffness matrix A and load vector b.

%Initialize A to have two more rows and columns.
A = zeros(M+2);
b = zeros(M,1);

for i = 1:M 
    if rem(i,2) ~= 0                   % if i is odd
        A(i,i) = A(i,i) + 1/(3*h);

        b(i) = b(i) + f(x(i)) * (-2)*h/3;

    elseif rem(i,2) == 0               % if i is even
        A(i,i) = A(i,i) + 3/(3*h);
        A(i+2,i) = A(i+2,i) - 8/(3*h);
        A(i,i+2) = A(i,i+2) - 8/(3*h);

        b(i) = b(i) + f(x(i)) * 8*h/3;
    end

    A(i+1,i) = A(i+1,i) + 2/(3*h);
    A(i,i+1) = A(i,i+1) - 3/(3*h);
end

%Remove two rows and columns form matrix A.
A = A(1:end-2, 1:end-2);

% Enforce the boundary conditions

A(1,1) = 1.e4;
A(end,end) = 1.e4;


function z = f(t)
% Stab f functionn - f(x) return x
z = t;
end

J'espère que le code se comportera comme prévu (la valeur que vous vous attendez à lire n'est pas si claire en dehors des limites de A)...

1
Rotem 9 févr. 2020 à 22:08