Comment changer le programme Matlab pour résoudre une équation avec la méthode des éléments finis?

Question

Je ne sais pas si cette question est plus liée aux mathématiques ou à la programmation et je suis absolument novice dans Matlab. Le programme FEM_50 applique la méthode des éléments finis à l'équation de Laplace -Uxx (x, y) - Uyy (x, y) = F (x, y) dans Omega. Comment le changer pour appliquer FEM à l'équation -Uxx (x, y) - Uyy (x, y) + U (x, y) = F (x, y)? Sur cette page: http://sc.fsu.edu/~burkardt/m_src/fem_50/fem_50.html fichiers de code supplémentaires au cas où vous en avez besoin.

function fem_50 () 

%% FEM_50 applies the finite element method to Laplace's equation. 
% 
% Discussion: 
% 
% FEM_50 is a set of MATLAB routines to apply the finite 
% element method to solving Laplace's equation in an arbitrary 
% region, using about 50 lines of MATLAB code. 
%  
% FEM_50 is partly a demonstration, to show how little it 
% takes to implement the finite element method (at least using 
% every possible MATLAB shortcut.) The user supplies datafiles 
% that specify the geometry of the region and its arrangement 
% into triangular and quadrilateral elements, and the location 
% and type of the boundary conditions, which can be any mixture 
% of Neumann and Dirichlet. 
%  
% The unknown state variable U(x,y) is assumed to satisfy 
% Laplace's equation: 
%  -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega 
% with Dirichlet boundary conditions 
%  U(x,y) = U_D(x,y) on Gamma_D 
% and Neumann boundary conditions on the outward normal derivative: 
%  Un(x,y) = G(x,y) on Gamma_N 
% If Gamma designates the boundary of the region Omega, 
% then we presume that 
%  Gamma = Gamma_D + Gamma_N 
% but the user is free to determine which boundary conditions to 
% apply. Note, however, that the problem will generally be singular 
% unless at least one Dirichlet boundary condition is specified. 
%  
% The code uses piecewise linear basis functions for triangular elements, 
% and piecewise isoparametric bilinear basis functions for quadrilateral 
% elements. 
%  
% The user is required to supply a number of data files and MATLAB 
% functions that specify the location of nodes, the grouping of nodes 
% into elements, the location and value of boundary conditions, and 
% the right hand side function in Laplace's equation. Note that the 
% fact that the geometry is completely up to the user means that 
% just about any two dimensional region can be handled, with arbitrary 
% shape, including holes and islands. 
%  
    clear 
% 
% Read the nodal coordinate data file. 
% 
    load coordinates.dat; 
% 
% Read the triangular element data file. 
% 
    load elements3.dat; 
% 
% Read the quadrilateral element data file. 
% 
    load elements4.dat; 
% 
% Read the Neumann boundary condition data file. 
% I THINK the purpose of the EVAL command is to create an empty NEUMANN array 
% if no Neumann file is found. 
% 
    eval ('load neumann.dat;', 'neumann=[];'); 
% 
% Read the Dirichlet boundary condition data file. 
% 
    load dirichlet.dat; 

    A = sparse (size(coordinates,1), size(coordinates,1)); 
    b = sparse (size(coordinates,1), 1); 
% 
% Assembly. 
% 
    for j = 1 : size(elements3,1) 
    A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ... 
     + stima3(coordinates(elements3(j,:),:)); 
    end 

    for j = 1 : size(elements4,1) 
    A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ... 
     + stima4(coordinates(elements4(j,:),:)); 
    end 
% 
% Volume Forces. 
% 
    for j = 1 : size(elements3,1) 
    b(elements3(j,:)) = b(elements3(j,:)) ... 
     + det([1,1,1; coordinates(elements3(j,:),:)']) * ... 
     f(sum(coordinates(elements3(j,:),:))/3)/6; 
    end 

    for j = 1 : size(elements4,1) 
    b(elements4(j,:)) = b(elements4(j,:)) ... 
     + det([1,1,1; coordinates(elements4(j,1:3),:)']) * ... 
     f(sum(coordinates(elements4(j,:),:))/4)/4; 
    end 
% 
% Neumann conditions. 
% 
    if (~isempty(neumann)) 
    for j = 1 : size(neumann,1) 
     b(neumann(j,:)) = b(neumann(j,:)) + ... 
     norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ... 
     g(sum(coordinates(neumann(j,:),:))/2)/2; 
    end 
    end 
% 
% Determine which nodes are associated with Dirichlet conditions. 
% Assign the corresponding entries of U, and adjust the right hand side. 
% 
    u = sparse (size(coordinates,1), 1); 
    BoundNodes = unique (dirichlet); 
    u(BoundNodes) = u_d (coordinates(BoundNodes,:)); 
    b = b - A * u; 
% 
% Compute the solution by solving A * U = B for the remaining unknown values of U. 
% 
    FreeNodes = setdiff (1:size(coordinates,1), BoundNodes); 

    u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes); 
% 
% Graphic representation. 
% 
    show (elements3, elements4, coordinates, full (u)); 

    return 
end 

Answer 1

Je ne vais pas faire vos devoirs, ... mais voici quelques conseils:

1. Retour aux formules et travailler la forme faible de la nouvelle équation à résoudre, vous Vous noterez une nouvelle contribution à la matrice A, que vous devez également prendre en compte dans la section "Assemblage".

2. Trouver la matrice locale (élément par élément) expression de cette nouvelle contribution, vous pouvez l'ajouter aussi bien dans la section « Assemblée ». Vous verrez que c'est en fait la matrice de masse de l'élément.

3. Modifier la section "Assemblée" pour tous les types d'éléments (triangles et Quads)

4. et le tour est joué ...

Ce fut la façon "propre" de le faire. Il y a en fait une autre manière qui n'exigerait presque aucune modification du programme fourni. Il suffit de le transformer en une fonction Matlab pour vous permettre de résoudre à plusieurs reprises le problème d'origine avec différentes forces du corps (F (x, y)).

Ensuite, appelez à plusieurs reprises cette fonction lorsque vous modifiez la force du corps à:

New_body_Force = Original_body_force - Solution_from_previous_call 

Cela devrait converger vers le résultat souhaité.

Cette deuxième option est aussi élégante que la façon « propre » d'abord présenté mais je l'amour itérations de point fixe.

Espérons que cela aide.