Finite Differenzenmethode

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

thalia27

Forum-Anfänger


	Beiträge: 11

	Anmeldedatum: 14.04.20

	Wohnort: ---

	Version: ---

Verfasst am: 09.01.2021, 19:47 Titel: Finite Differenzenmethode

Hallo zusammen,

Ich habe versucht die finite Differenzenmethode in MATLAB zu lösen. Dabei ist mir aufgefallen, dass die Randbereiche nicht korrekt im Einheitsquadrat (0<x<1, 0<y<1) dargestellt werden und das die Funktion F(x,y) sich über die Auflösungsgröße verändert, was ja eigentlich nicht sein sollte.

Code:

clc
clear

n = 50; %Anzahl der Bereiche
h = 1/n; %Schrittweite ein x- und y-Richtung
A = Koeffizientenmatrix(n);

RB = Randbedingungen(n,h);
%Ermitteln des Vektors Q
Q=VektorQ(n);

%Berechnen des Vektors b aus Q+RB
b=Q+RB;

%Ermitteln des Vektors U zur spÃ¤teren Darstellung
U=VektorU(A, b, n);

%Ermitteln von U durch Gauss- Elimination
UGauss = Gauss(A,b);

%Zusammenstellen von RB und U
F = MatrixF (U, n, h);

%% Koeffizientenmatrix

function [Koeffizientenmatrix]=Koeffizientenmatrix(n)

I = eye(n-1)*(-1); %negierte Einheitsmatrix erstellen

B = zeros(n-1); %Nullmatrix B erstellen
for i = 1: n-1 %Die zu betrachtende Punkte mit einer 4 ausfÃ¼llen
B(i,i) = 4;
end

A = zeros((n-1)^2); %Nullmatrix A mit ((n-1)*(n-1))^2 EintrÃ¤gen erstellen
for i=0 : (n-2) %Die Diagonale von A mit B ausfÃ¼llen
A(i*(n-1)+1:(i+1)*(n-1),i*(n-1)+1:(i+1)*(n-1))=B;
end

for i=0 : (n-2)-1 %Rechts neben der Diagonalmatrix mit I (Einheitsmatrix) ausfÃ¼llen
A(i*(n-1)+1:(i+1)*(n-1),(i+1)*(n-1)+1:(i+2)*(n-1))=I;
end

for i=0 : (n-2)-1 %Links neben der Diagonalmatrix mit I (Einheitsmatrix) ausfÃ¼llen
A((i+1)*(n-1)+1:(i+2)*(n-1),i*(n-1)+1:(i+1)*(n-1))=I;
end

for i=1 : ((n-1)^2)-1 %-1 links und rechts der Hauptdiagonalen eintragen
A(i,i+1) = -1;
A(i+1,i) = -1;
end
Koeffizientenmatrix = A;

end
%% Randbedingungen

function [Randbedingungen]= Randbedingungen(n,h)

%Randbedingungen der 4 Randbereiche formulieren
RBu = zeros((n-1)^2,1);
RBo = zeros((n-1)^2,1);
RBli = zeros((n-1)^2,1);
RBre = zeros((n-1)^2,1);
% RB fÃ¼r die einzelnen definieren

%unten
for i = 1 : n-1
RBu(i) = i*h;
end

%oben
for i = 1 : n-1
RBo(i+(n-1)*(n-2)) = 2*RBu(i);
end

%linke RB = 0, keine Elemente einzufÃ¼gen

%rechts
for i = 1 : n-1
RBre(i*(n-1)) = 1 + i*h;
end

%Matrizen fÃ¼r RB zusammenfÃ¼hren
RB = RBu + RBo+ RBli + RBre;

Randbedingungen = RB;

end

%%%%%% Zuordnung eines Vektors in eine Matrix%%%%%%
function F = MatrixF (U, n, h)
tic %Zeit messen
%idxmat = reshape(uint64(U), n-1, n-1); %geringere SpeicherkapazitÃ¤t
idxmat = reshape(U, n-1, n-1); %hÃ¶here SpeicherkapazitÃ¤t
Randli = [1:-h:0]'; %RB links in Matrix einbauen
Randre = [0:h:1]'; %RB rechts in Matrix einbauen
Randun = [0+h:h:1-h];%RB unten in Matrix einbauen
Randob = [1-h:-h:0+h];%RB oben in Matrix einbauen
T = [Randob;idxmat;Randun];%Zusammenstellen von Matrix U und RBoben RBunten
F = [Randli T Randre];%Zusammenstellen aller RB und U
xyaxis= [0:h:1];%Skalierung
contour3 (xyaxis, xyaxis, F);%Grafische Darstellung
toc %Zeit messen
save('Probe.mat') %Konventieren von .m Datei in .mat %Dateiname muss angepasst werden
S = whos('-file','Probe.mat'); %Dateiname muss angepasst werden
for k = 1:length(S)
disp(['' S(k).name ...
'' mat2str(S(k).size) ...
'' S(k).class]);
end
end

function [VektorU] = VektorU(A, b, n)

%% Funktionsparameter definieren
maxit = 1000;
x = zeros((n-1)*(n-1),1);
%x = rand((n-1)*(n-1),1);
tol = 1e-5;

%% Erstellung der x-Werte (Jacobi-Verfahren)
%dx=zeros((n-1)*(n-1),1);
for i=1:maxit
summe=0;
for k=1:(n-1)*(n-1)
dx(k)=b(k);
for j=1:(n-1)*(n-1)
dx(k)=dx(k)-A(k,j)*x(j);
end
dx(k)=dx(k)/A(k,k);
x(k)=x(k)+dx(k);
if (dx(k) <= 0)
summe = summe - dx(k);
else
summe = summe + dx(k);
end
end
if(summe <= tol)
break
end
end
%% Ausgabe des Vektors x als Funktionswert VektorU
VektorU=x;
end

function [VektorQ]= VektorQ(n)

%% Erzeugen von x und y Werten
[x,y]= meshgrid(0 :(1/n): 1,1:-(1/n):0);

%% Ergebnismatrix fÃ¼r Piosson Gleichung mit x- und y-Werten
F = zeros(n-1,n-1);

for j= 1 : (n-1)
for k= 1:(n-1)
F(j,k)= sin(4*pi*x(j+1,k+1))*cos(4*pi* y(j+1,k+1));
end
end

%% Vektor Q aus Matrix F erzeugen
Q = zeros((n-1)*(n-1),1);
m = 1;

for j = 1 : (n-1)
for k= 1:(n-1)
Q(m,1)= F(j,k);
m=m+1;
end
end

%% Einpflegen der Randbediungen in in Vektor Q + Ausgabe
VektorQ=Q;

end

function [Gauss] = Gauss(A,b)
UGauss = A\b;
Gauss = UGauss;
end

Funktion ohne Link?

ich denke, dass das Bild so aussehen sollte. Stattdessen liefert das Programm dieses Bild

Kann jemand mir sagen, wo der Fehler ist?
Vielen Dank


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.