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

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

  
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Gleichungssystem lösen
 

Henne24

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 09.11.2011, 14:25     Titel: Gleichungssystem lösen
  Antworten mit Zitat      
Hallo, ich habe ein Problem beim definieren eines GLS und das sich anschließende Lösen mit fsolve.

Hier mal mein Quellcode:
Code:
%Berechnung der Füllstände im Behälterverbund BH51 und BH52
%Konstanten
D=12.0/1000.0; % Rohrdurchmesser in m
k=0.03/1000.0; % relative Rohrrauhigkeit in m
zetaKr=1.33;   % Zetawert für 90°-Rohrkrümmer
g=9.81;        % Erdbeschleunigung
rho=1000.0;    % Dichte
PI=3.14;       % PI
%Parameter
LaengeR1=0.3; %m
LaengeR2=1.7; %m
AnzKR1=2.0;
AnzKR2=5.0;
Grundfl=0.25*0.25;
dh52max=0.55 + 0.17;
FuellhoeheB52=0.1;%m
FuellhoeheB51=0.17;%m

lambda=sqrt(1.0/(1.14-2.0*log10(k/D))); % Durchflussbeiwert

%GLS

h52n=((FuellhoeheB52*Grundfl)-Q1+Q2)/Grundfl;
pt52=rho*g*h52n;
pvR52=(8*lambda*LaengeR1*rho*Q1*abs(Q1))/D^5*PI^2*g;
pvE52=(8*rho*Q1*abs(Q1)*zetaKr*AnzKR1)/D^4*PI^2*g;
0=-pt52+(pvR52+pvE52;

h51n=((FuellhoeheB51*Grundfl)+Q1-Q2)/Grundfl;
pt51=-2E10*Q2*Q2+833333.0*Q2+25000.0;
dh51=rho*g*(dh52max-h51n);
pvR51=(8*lambda*LaengeR1*rho*Q2*abs(Q2))/D^5*PI^2*g;
pvE51=(8*rho*Q1*abs(Q1)*zetaKr*AnzKR2)/D^4*PI^2*g;
0=-pt51+(dh51+pvR51+pvE52;
 


Das Gleichungssystem soll nach Q1 und Q2 gelöst werden. Also es soll die Nullstelle gefunden werden. Vielleicht kann man ja auch mit
Code:
eine Funktion definieren. Ich habe aber keine Ahnung mehr wie man das programmiert. Das ist schon ne Weile her.
Ich würde mich sehr freuen wenn Ihr mir helfen könntet.

Stefan


Harald
Forum-Meister

Forum-Meister

Beiträge: 24.502
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 09.11.2011, 15:52     Titel:
  Antworten mit Zitat      
Hallo,

ich würde dafür FSOLVE verwenden.

Code:
[loesung, abweichung] = fsolve(@gleichung, startpunkt)
wobei gleichung.m wie folgt aussieht:

Code:
function y = gleichung(x)
Q1 = x(1);
Q2 = x(2);
%Berechnung der Füllstände im Behälterverbund BH51 und BH52
%Konstanten
D=12.0/1000.0; % Rohrdurchmesser in m
k=0.03/1000.0; % relative Rohrrauhigkeit in m
zetaKr=1.33;   % Zetawert für 90°-Rohrkrümmer
g=9.81;        % Erdbeschleunigung
rho=1000.0;    % Dichte
PI=3.14;       % PI
%Parameter
LaengeR1=0.3; %m
LaengeR2=1.7; %m
AnzKR1=2.0;
AnzKR2=5.0;
Grundfl=0.25*0.25;
dh52max=0.55 + 0.17;
FuellhoeheB52=0.1;%m
FuellhoeheB51=0.17;%m

lambda=sqrt(1.0/(1.14-2.0*log10(k/D))); % Durchflussbeiwert

%GLS

h52n=((FuellhoeheB52*Grundfl)-Q1+Q2)/Grundfl;
pt52=rho*g*h52n;
pvR52=(8*lambda*LaengeR1*rho*Q1*abs(Q1))/D^5*PI^2*g;
pvE52=(8*rho*Q1*abs(Q1)*zetaKr*AnzKR1)/D^4*PI^2*g;
y(1) =-pt52+(pvR52+pvE52;

h51n=((FuellhoeheB51*Grundfl)+Q1-Q2)/Grundfl;
pt51=-2E10*Q2*Q2+833333.0*Q2+25000.0;
dh51=rho*g*(dh52max-h51n);
pvR51=(8*lambda*LaengeR1*rho*Q2*abs(Q2))/D^5*PI^2*g;
pvE51=(8*rho*Q1*abs(Q1)*zetaKr*AnzKR2)/D^4*PI^2*g;
y(2) =-pt51+(dh51+pvR51+pvE52;


Bitte überprüfe die Klammern im Code. Zudem erscheint es mir merkwürdig, dass du pvE51 am Ende definierst, aber nicht verwendest.

Grüße,
Harald
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.502
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 09.11.2011, 15:58     Titel:
  Antworten mit Zitat      
P.S.:
1. pi ist in MATLAB vordefiniert.
2. Klammernsetzung z.B. bei
pvR52=(8*lambda*LaengeR1*rho*Q1*abs(Q1))/D^5*PI^2*g;
Willst du wirklich nur durch D^5 teilen, nicht durch D^5*PI^2*g? Wenn du durch letzteres teilen willst, brauchst du Klammern um den Nenner. Bei den anderen Ausdrücken natürlich entsprechend.

P.P.S.:
Da die Gleichung um Q1=Q2=0 nicht differenzierbar ist, liefern fminsearch evtl. bessere Ergebnisse:
Code:
[loesung, abweichung] = fminsearch(@(x) norm(gleichung123(x)), startwert)
Private Nachricht senden Benutzer-Profile anzeigen
 
Henne24

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 10.11.2011, 12:29     Titel:
  Antworten mit Zitat      
Vielen Dank für die Hinweise. Ich hab es jetzt so gelöst.
Code:
function y = MaGl(x)

%Konstanten
D=12.0/1000.0; % Rohrdurchmesser in m
k=0.03/1000.0; % relative Rohrrauhigkeit in m
zetaKr=1.33;   % Zetawert für 90°-Rohrkrümmer
g=9.81;        % Erdbeschleunigung
rho=1000.0;    % Dichte

%Parameter
LaengeR1=0.3; %m
LaengeR2=1.7; %m
AnzKR1=2.0;
AnzKR2=5.0;
Grundfl=0.25*0.25;
dh52max=0.55 + 0.17;
FuellhoeheB52=0.1;%m
FuellhoeheB51=0.17;%m

lambda=sqrt(1.0/(1.14-2.0*log10(k/D))); % Durchflussbeiwert
   
h52n=((FuellhoeheB52*Grundfl)-x(1)+x(2))/Grundfl;
pt52=rho*g*h52n;
pvR52=(8*lambda*LaengeR1*rho*x(1)*abs(x(1)))/(D^5*pi^2*g);
pvE52=(8*rho*x(1)*abs(x(1))*zetaKr*AnzKR1)/(D^4*pi^2*g);

h51n=((FuellhoeheB51*Grundfl)+x(1)-x(2))/Grundfl;
pt51=-2E10*x(2)*x(2)+833333.0*x(2)+25000.0;
dh51=rho*g*(dh52max-h51n);
pvR51=(8*lambda*LaengeR2*rho*x(2)*abs(x(2)))/(D^5*pi^2*g);
pvE51=(8*rho*x(1)*abs(x(1))*zetaKr*AnzKR2)/(D^4*pi^2*g);

y(1)=-pt52+(pvR52+pvE52);
y(2)=-pt51+(dh51+pvR51+pvE51);


Und dann aus einem anderen m-file mit fsolve die Lösung gefunden.
Code:
format long
options(2)=1e-12;
z=fsolve('MaGl',[1;1],options);


Stefan
 
Henne24

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 15.11.2011, 11:36     Titel:
  Antworten mit Zitat      
Hallo, ich möchte die Funktion jetzt so anpassen, dass die Fuellhoehen auch von außen an die Funktion übergeben werden können und ich das Gleichungssystem so wie bisher mit fsolve lösen kann.
Ich habe schon diese Möglichkeit probiert.
Code:
   function y = Stroeme(x, FuellhoeheB51,FuellhoeheB52)
   
    %Parameter
    LaengeR1=0.3; %m
    LaengeR2=1.7; %m
    AnzKR1=2.0;
    AnzKR2=5.0;
    dh52max=0.69;

    %Konstanten
    D=12.0/1000.0; % Rohrdurchmesser in m
    k=0.03/1000.0; % relative Rohrrauhigkeit in m
    zetaKr=1.33;   % Zetawert für 90°-Rohrkrümmer
    g=9.81;        % Erdbeschleunigung
    rho=1000.0;    % Dichte
   
    lambda=sqrt(1.0/(1.14-2.0*log10(k/D))); % Durchflussbeiwert    
    pt52=rho*g*FuellhoeheB52;
    pvR52=(8*lambda*LaengeR1*rho*x(1)*abs(x(1)))/(D^5*pi^2*g);
    pvE52=(8*rho*x(1)*abs(x(1))*zetaKr*AnzKR1)/(D^4*pi^2*g);

    pt51=-2E10*x(2)*x(2)+833333.0*x(2)+25000.0;
    ptdh51=rho*g*(dh52max-FuellhoeheB51);
    pvR51=(8*lambda*LaengeR2*rho*x(2)*abs(x(2)))/(D^5*pi^2*g);
    pvE51=(8*rho*x(2)*abs(x(2))*zetaKr*AnzKR2)/(D^4*pi^2*g);
   
    y(1)=-pt52+(pvR52+pvE52);
    y(2)=-pt51+(ptdh51+pvR51+pvE51);
    end

Wie kann ich die Funktion dann mit fsolve lösen und die Fuellhoehen an die Funktion Stroeme ubergeben?

Code:
z=fsolve('Stroeme',[1;1],options);


Stefan
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.502
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 15.11.2011, 11:42     Titel:
  Antworten mit Zitat      
Hallo,

über anonymous Function Handles (evtl. auch in der Hilfe dazu nachlesen)

Code:
FuellhoeheB51 = ...
FuellhoeheB52 = ...
z=fsolve(@(x) Stroeme(x, FuellhoeheB51, FuellhoeheB52) ,[1;1],options);  


Grüße,
Harald

P.S.: wenn zu einem "beantworteten" Thread weitere Fragen bestehen, bitte den Status auf "offen" setzen.
Private Nachricht senden Benutzer-Profile anzeigen
 
Henne24

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 15.11.2011, 12:05     Titel:
  Antworten mit Zitat      
Ich kann mich nur wiederum bei dir bedanken. Es funktioniert. Hier nochmal meine Lösung:
Code:
   function y = Stroeme(x,FuellhoeheB51,FuellhoeheB52)
   
    %Parameter
    LaengeR1=0.3; %m
    LaengeR2=1.7; %m
    AnzKR1=2.0;
    AnzKR2=5.0;
    dh52max=0.69;
    %Konstanten
    D=12.0/1000.0; % Rohrdurchmesser in m
    k=0.03/1000.0; % relative Rohrrauhigkeit in m
    zetaKr=1.33;   % Zetawert für 90°-Rohrkrümmer
    g=9.81;        % Erdbeschleunigung
    rho=1000.0;    % Dichte
   
    lambda=sqrt(1.0/(1.14-2.0*log10(k/D))); % Durchflussbeiwert    
    pt52=rho*g*FuellhoeheB52
    pvR52=(8*lambda*LaengeR1*rho*x(1)*abs(x(1)))/(D^5*pi^2*g);
    pvE52=(8*rho*x(1)*abs(x(1))*zetaKr*AnzKR1)/(D^4*pi^2*g);

    pt51=-2E10*x(2)*x(2)+833333.0*x(2)+25000.0;
    ptdh51=rho*g*(dh52max-FuellhoeheB51);
    pvR51=(8*lambda*LaengeR2*rho*x(2)*abs(x(2)))/(D^5*pi^2*g);
    pvE51=(8*rho*x(2)*abs(x(2))*zetaKr*AnzKR2)/(D^4*pi^2*g);
   
    y(1)=-pt52+(pvR52+pvE52);
    y(2)=-pt51+(ptdh51+pvR51+pvE51);
    end


Und hier der Aufruf des Solvers:
Code:
function Z = gleichung5152(u1,u2)
FuellhoeheB51=u1;
FuellhoeheB52=u2;
Z = fsolve(@(x)Stroeme(x,FuellhoeheB51,FuellhoeheB52),[1,1]);
end
 


Das hätte ich nie herausgefunden. Es wäre ja schon eine große Hilfe wenn die Hilfe auf deutsch wäre.

Also vielen Dank nochmal
 
Neues Thema eröffnen Neue Antwort erstellen



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 | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2025 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.