Matroids Matheplanet Forum Index
Moderiert von mire2
Mathematische Software & Apps » Matlab » Adams-Moulton-Verfahren
Autor
Universität/Hochschule Adams-Moulton-Verfahren
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Themenstart: 2021-11-30

https://www.matheplanet.com/matheplanet/nuke/html/uploads/b/54010_Unbenltkgjt_lgtannt.PNG Hallo zusammen, könnte mir jemand bitte dabei helfen? Ich muss diese Formel auf Matlab implementieren : https://www.matheplanet.com/matheplanet/nuke/html/uploads/b/54010_Unbentogijtoigjtoannt.PNG Und das ist die Tabelle für beta und k = 1,2,3 : https://www.matheplanet.com/matheplanet/nuke/html/uploads/b/54010_Unbekhglgtjltgtnannt.PNG Mein Problem ist dass das Verfahren implizit ist und weiß nicht was ich statt f(x_1,y_1) einsetzen soll Danke im Voraus!


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.1, vom Themenstarter, eingetragen 2021-11-30

https://www.matheplanet.com/matheplanet/nuke/html/uploads/b/54010_71472e10-1008-477c-b9e6-abd4f7da3e1f.jpg


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4607
  Beitrag No.2, eingetragen 2021-11-30

\quoteon(2021-11-30 20:32 - s-amalgh im Themenstart) Mein Problem ist dass das Verfahren implizit ist und weiß nicht was ich statt f(x_1,y_1) einsetzen soll \quoteoff Du musst ausnutzen, dass hier $f(x_i,y_i)=A\,y_i$ linear in $y_i$ ist. Wenn du das in die Gleichung, die $y_{i+1}$ implizit definiert, einsetzt, erhältst du ein lineares Gleichungssystem für $y_{i+1}$. --zippy


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.3, vom Themenstarter, eingetragen 2021-11-30

Danke erstmal für deine Antwort! Wie kann ich das in Matlab implementieren ? Könntest du mir bitte bei dem ersten Schritt helfen? Dnke im Voraus!


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4607
  Beitrag No.4, eingetragen 2021-11-30

\quoteon(2021-11-30 22:37 - s-amalgh in Beitrag No. 3) Wie kann ich das in Matlab implementieren ? \quoteoff Matlab kann lineare Gleichungssysteme lösen, das hast du schonmal verwendet: \quoteon(2021-11-06 16:25 - s-amalgh in Beitrag No. 24) \sourceon Matlab z = -m\z \sourceoff \quoteoff Hier lautet das Gleichungssystem $\left(1-\frac h2A\right)y_1 = \left(1+\frac h2A\right)y_0$, wobei $1$ die $2\times2$-Einheitsmatrix bezeichnet.


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.5, vom Themenstarter, eingetragen 2021-11-30

wie kommst du auf (1-h/2*A) ?


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4607
  Beitrag No.6, eingetragen 2021-11-30

\quoteon(2021-11-30 22:50 - s-amalgh in Beitrag No. 5) wie kommst du auf (1-h/2*A) ? \quoteoff Die ursrpüngliche Gleichung lautet wegen $\beta_0^*=\beta_1^*=\frac12$$$ y_1 = y_0 + \frac h2\Bigl[f(x_1,y_1)+f(x_0,y_0)\Bigr] = y_0 + \frac h2\Bigl[A\,y_1+A\,y_0\Bigr] \;, $$und Sortieren nach $y_1$ und $y_0$ liefert$$ y_1-\frac h2A\,y_1 = y_0+\frac h2A\,y_0 \quad\iff\quad \left(1-\frac h2A\right)y_1 = \left(1+\frac h2A\right)y_0 \;. $$


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.7, vom Themenstarter, eingetragen 2021-11-30

Achso ok.. könntest du mir bitte zeigen wie ich das in matlab implementieren kann? nur den ersten Schritt, ich ach dann der Rest allein


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4607
  Beitrag No.8, eingetragen 2021-11-30

\quoteon(2021-11-30 22:37 - s-amalgh in Beitrag No. 3) Wie kann ich das in Matlab implementieren ? \quoteoff Unter der Annahme, dass a die Matrix $A$, y den Vektor $y_0$ und h die Schrittweite $h$ enthält, ergibt sich $y_1$ als: \sourceon matlab (eye(2)-.5*h*a)\((eye(2)+.5*h*a)*y) \sourceoff


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.9, vom Themenstarter, eingetragen 2021-12-01

\sourceon Matlab function t1 y0 = [1; 0]; a = 0; b = 10; h = [1e-1,1e-2,1e-3]; A = [0, 1; -4*pi^2, 0]; f = @(x,y) (A*y); for k = 1:3 for j = 1:length(h) [~,y] = adammoult(a, b,y0,h(j),1,f); end end function [x,y] = adammoult(a,b,y0,h,k,rhs) beta = [1/2 1/2 0 0; ... % k = 1 5/12 8/12 -1/12 0; ... % k = 2 9/24 19/24 -5/24 1/24]; ... % k = 3 x = a:h:b; N = length(x); n = length(y0); y = [y0(:), zeros(n,N-1)]; if k == 1 for i = 1:N-1 y(:,i+1) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,i)); end elseif k ==2 for i = 2:N-1 y(:,2) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,1)); % Hier überlege ich noch was ich schreiben soll % end end end \sourceoff So ? noch nicht fertig aber als idee


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4607
  Beitrag No.10, eingetragen 2021-12-01

Eine Fallunterscheidung nach dem Wert von $k$ brauchst du eigentlich nicht. Du kannst für alle $k$: 1. Die linke Seite von ...\... aus $\beta_0^*$ sofort berechnen. 2. Die rechte Seite von ...\... in einer Schleife als Summe berechnen.


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.11, vom Themenstarter, eingetragen 2021-12-01

wie meinst du das ? Mein bisheriger Code ist so : \sourceon Matlab function t1 y0 = [1; 0]; a = 0; b = 10; h = [1e-1,1e-2,1e-3]; A = [0, 1; -4*pi^2, 0]; for k = 1:3 for j = 1:length(h) [~,y] = adammoult(a, b,y0,h(j),k,A); end end function [x,y] = adammoult(a,b,y0,h,k,A) beta = [1/2 1/2 0 0; ... % k = 1 5/12 8/12 -1/12 0; ... % k = 2 9/24 19/24 -5/24 1/24]; ... % k = 3 x = a:h:b; N = length(x); n = length(y0); y = [y0(:), zeros(n,N-1)]; if k == 1 for i = 1:N-1 y(:,i+1) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,i)); end elseif k == 2 y(:,2) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,1)); for i = 2:N-1 y(:,i+1) = (eye(2)-5/12*h*A)\((eye(2)+8/12*h*A)*y(:,i) - 1/12*h*A*y(:,i-1)); end elseif k == 3 y(:,2) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,1)); y(:,3) = (eye(2)-5/12*h*A)\((eye(2)+8/12*h*A)*y(:,2) - 1/12*h*A*y(:,1)); for i = 3:N-1 y(:,i+1) = (eye(2)-9/24*h*A)\((eye(2)+19/24*h*A)*y(:,i) - 5/24*h*A*y(:,i-1) + 1/24*h*A*y(:,i-2)); end end end end \sourceoff


   Profil
zippy
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 24.10.2018
Mitteilungen: 4607
  Beitrag No.12, eingetragen 2021-12-01

Das sieht aus, als könnte es funktionieren (ich habe es mit aber nicht im Detail angeschaut). Allerdings programmierst du einiges aus, was du durch Schleifen oder Funktionen ersetzen könntest. Ein Symptom dafür ist, dass du die Matrix beta gar nicht verwendest.


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.13, vom Themenstarter, eingetragen 2021-12-01

Könntst du bitte mein Code bearbeiten?


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.14, vom Themenstarter, eingetragen 2021-12-01

wie kann ich jetzt plotten? plot(h,y(1,:)) so? wenn ich run mache zeigt dass die beide Vektoren die gleiche Länge sein müssen. Wie mach ich das?


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.15, vom Themenstarter, eingetragen 2021-12-01

Hast du idee?


   Profil
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 11600
Wohnort: Wien
  Beitrag No.16, eingetragen 2021-12-01

Hallo s-amalgh, überlege Dir, welche Bedeutung $h$ und der derzeit nicht verwendete Vektor x haben. Servus, Roland


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.17, vom Themenstarter, eingetragen 2021-12-01

Ich habe nicht verstanden was du gemeint hast.. Kannst du mir bitte das mehr erklären?


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.18, vom Themenstarter, eingetragen 2021-12-01

Meinst du das so? für k = 1: \sourceon Matlab h = [1e-1,1e-2,1e-3]; #for k = 1:3 for j = 1:length(h) [~,y] = adammoult(a, b,y0,h(j),1,A) plot(0:h(j):10,y(1,:)) hold on end \sourceoff


   Profil
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 11600
Wohnort: Wien
  Beitrag No.19, eingetragen 2021-12-01

Hallo s-amalgh, die Funktion plot() erwartet zwei Vektoren mit den x- und y-Koordinaten der zu zeichnenden Punkte. Die Zahl $h$ ist die Schrittweite, also der horizontale Abstand zweier benachbarter Punkte. Die Funktion adammoult() liefert zwei Ergebnisse zurück, von denen Du das erste nicht verwendest. Ich hoffe, das hilft Dir, Roland [Die Antwort wurde nach Beitrag No.17 begonnen.]


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.20, vom Themenstarter, eingetragen 2021-12-01

für k = 1: \sourceon Matlab h = [1e-1,1e-2,1e-3]; #for k = 1:3 for j = 1:length(h) [~,y] = adammoult(a, b,y0,h(j),1,A) plot(0:h(j):10,y(1,:)) hold on end \sourceoff Meinst du das so oder wie?


   Profil
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 11600
Wohnort: Wien
  Beitrag No.21, eingetragen 2021-12-01

Hallo s-amalgh, nein, ich meinte, dass Du den Vektor x, den Du hoffentlich nicht ohne Grund in der Funktion adammoult berechnen lässt, verwendest wofür er gedacht ist. Servus, Roland


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.22, vom Themenstarter, eingetragen 2021-12-01

Du meintest das so oder: Für k = 1: \sourceon Matlab h = [1e-1,1e-2,1e-3]; for j = 1:length(h) [x,y] = adammoult(a, b,y0,h(j),1,A) plot(x,y(1,:)) hold on end \sourceoff


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.23, vom Themenstarter, eingetragen 2021-12-01

Mein Code : \sourceon Matlab function Alghabra07 function [x,y] = adammoult(a,b,y0,h,k,A) % Gitterpunkte der N"aherungen x = a:h:b; % Anzahl der Gitterpunkte N = length(x); % Dimension der Dgl n = length(y0); % Initialisierung der N"aherungen y = [y0(:), zeros(n,N-1)]; %Adams-Moulton-Verfahren % K = 1 if k == 1 for i = 1:N-1 y(:,i+1) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,i)); end % K = 2 elseif k == 2 y(:,2) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,1)); for i = 2:N-1 y(:,i+1) = (eye(2)-5/12*h*A)\((eye(2)+8/12*h*A)*y(:,i) - 1/12*h*A*y(:,i-1)); end % K = 3 elseif k == 3 y(:,2) = (eye(2)-.5*h*A)\((eye(2)+.5*h*A)*y(:,1)); y(:,3) = (eye(2)-5/12*h*A)\((eye(2)+8/12*h*A)*y(:,2) - 1/12*h*A*y(:,1)); for i = 3:N-1 y(:,i+1) = (eye(2)-9/24*h*A)\((eye(2)+19/24*h*A)*y(:,i) - 5/24*h*A*y(:,i-1) + 1/24*h*A*y(:,i-2)); end end end % Beispiel y0 = [1; 0]; % x_0 a = 0; % x_end b = 10; % Schrittweiten h = [1e-1,1e-2,1e-3]; % Matrix A = [0, 1; -4*pi^2, 0]; for k = 1:3 for j = 1:length(h) [x,y] = adammoult(a, b,y0,h(j),k,A) figure(k); plot(x,y(1,:)); set(gca,'FontSize',24); title(['k = ',num2str(k)]); hold on; end hold off; legend({'h(1)','h(2)','h(3)'},'Location','southwest'); legend('boxoff'); end end \sourceoff Und das sind die Plots : https://www.matheplanet.com/matheplanet/nuke/html/uploads/b/54010_Unbenkflghlbvghlitbtgrtannt.PNG Habe ich irgendein Fehler?


   Profil
rlk
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 16.03.2007
Mitteilungen: 11600
Wohnort: Wien
  Beitrag No.24, eingetragen 2021-12-05

Hallo s-amalgh, ich habe mir Dein Programm nicht im Detail angesehen, aber die Plots sehen so aus als wenn es funktionieren würde. In der Legende würde ich die Werte der verwendeten Schrittweite statt der wenig informativen Zeichenketten "h(1)" usw. ausgeben. Servus, Roland


   Profil
s-amalgh
Wenig Aktiv Letzter Besuch: vor mehr als 3 Monaten
Dabei seit: 16.12.2020
Mitteilungen: 374
  Beitrag No.25, vom Themenstarter, eingetragen 2021-12-05

Danke :) Könntest du mir bitte dabei helfen? https://www.matheplanet.com/matheplanet/nuke/html/viewtopic.php?topic=256608&start=0&lps=1863759#v1863759 Es wäre sehr nett von dir, da ich das heute abgeben muss.. Danke im Voraus!


   Profil
s-amalgh hat die Antworten auf ihre/seine Frage gesehen.
s-amalgh wird per Mail über neue Antworten informiert.

Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2023 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]