Differentialgleichungen < Matlab < Mathe-Software < Mathe < Vorhilfe
|
| Aufgabe | [mm] \phi ''=(4/(f^2*B^2))*((v/r*cos(\beta) [/mm] - [mm] a*sin(\beta))*(KG-T/2)*cos(\phi) [/mm] - [mm] g*GZ(\phi));
[/mm]
soll mit MatLab gelöst werden. Ziel ist, [mm] \phi [/mm] gegen die Zeit (t=0:1:500s) zu plotten
Anfangsbedingungen sind t=0; [mm] \phi=0; \phi [/mm] '=0 |
Ich habe diese Frage in keinem Forum auf anderen Internetseiten gestellt.
Hallo MatLab-Kenner!
Ich entwickle (manche erinnern sich vll an mein letztes Problem) immer noch ein Programm, mit dem ich das Krängungsverhalten eines Schiffes im Vollkreis untersuchen soll. Die Unterprogramme habe ich z.T. Dank dieses Forums hingekriegt und nun geht es also darum, diese m-files mit einem Hauptprogramm zusammenzufassen und auf diese Art und Weise eine Differentialgleichung zweiter Ordnung zu lösen.
g, f, B, KG und T sind feste Werte, die ich in einem Unterprogramm ('Schiff') zusammengefasst habe,
v, r, [mm] \beta [/mm] und a sind zeitabhängige Variable, die ich jede für sich in einem Unterprogramm mit einer spline-Approximation aus Stützstellen errechnet habe.
[mm] GZ(\phi)) [/mm] ist ebenfalls in einem Unterprogramm als function-file geschrieben (damit hab ich beim Zugreifen darauf auch noch Schwierigkeiten).
Auf all diese Unterprogramme (6 an der Zahl) soll nun das Hauptprogramm zugreifen und für die jeweiligen Zeitschritte die Werte in die DGL einsetzen.
Ich habe ziemlich wenig, um nicht zu sagen, gar keine, Erfahrung mit Schleifen... Könnt ihr mir Tipps geben? Ich habe erstmal ein Flussdiagramm erstellt, aber die Umsetzung gelingt mir einfach nicht. Mein Dozent sprach noch von der Taylor-Reihe, um uns einen Hinweis für die Bearbeitung zu geben... Tja... stimmt, sowas hatte ich in Mathe mal :-/ Aber wie ich das in MatLab schreibe...
Ich bin für jede Anregung dankbar!
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 19:07 Mi 16.05.2007 | | Autor: | cleophus |
Hallo nochmal!
Also ich war unterdesse natürlich nicht untätig und hab mich u.a. in diesem Forum über Schleifen schlau gelesen. Jetzt hab ich da eine zusammengebastelt, die nach einigen schlaflosen Nächten auch keine offensichtlichen Fehlermeldungen mehr produziert, aber ich habe ein großes Problem:
Ich will eigentlich n=500 berechnen lassen, und das errechnete [mm] \phi [/mm] dann plotten lassen. Allerdings sieht der Graph nur so lange ansatzweise brauchbar aus, wie ich n nicht größer werden lasse als ca 130. Kann mir das jemand erklären? Wie kann ich n größer machen?
% Unterprogramme einlesen
Geschwindigkeit;
Beschleunigung;
Kruemmung;
Driftwinkel;
% Konstanten
g=9.81;
[mm] c=4/(f^2*B^2);
[/mm]
Schiff;
% Anfangsbedingungen
dt=1; % Schrittweite in der Zeit
n=100; % Anzahl der Zeitschritte
phi=zeros(1,501);
phi1=zeros(1,501);
phi2=[];
% Algorithmus
% Schleife
for k=0:dt:n-1;
phi2(k+1)=c*((v(k+1).^2./r(k+1).*cos(beta(k+1))-a(k+1).*sin(beta(k+1)))*(KG-T/2)*cos(phi(k+1)) - g*GZ(phi(k+1)));
[mm] phi((k+2))=phi(k+1)+dt*phi1(k+1)+(dt^2)/2*phi2(k+1);
[/mm]
phi1((k+2))=phi1(k+1)+dt*phi2(k+1);
end
p=(pi/180).*phi
plot(t,p)
xlabel('Zeit t [s]')
ylabel('Krängungswinkel [mm] \phi [/mm] [°]')
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 17:46 Do 17.05.2007 | | Autor: | BKM |
Hallo
Könnten 'wir' vielleicht das vollständige .m File haben? Erspart u.a. das erneute eintippen ( jedenfalls was mich betrifft).
Beste Grüße
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 10:50 Sa 19.05.2007 | | Autor: | cleophus |
Ups, sorry, ich hatte die Hoffnung schon fast aufgegeben. Ich konnte aber einen Aufschub aushandeln, neues Fälligkeitsdatum ist nun erst Freitag, 25.05!
Wie gesagt sind es 6 Unterprogramme, eins davon ein function files. Ich gebe sie hier mal der Vollständigkeit halber alle an (vll kann ja jemand von euch die später noch anderweitig gebrauchen 
Schiff.m
%Schiffsdaten
global B T f KG KM
B=[21.5];
T=[7.96];
f=[0.8];
KG=[10.0];
KM=[11.0];
Beschleunigung.m
%Beschleunigung für MS Beluga Revolution
%loaded condition, deep water
%
global a
x=[0,115,244,355,482,500];
y=[6.173,3.7554,3.1381,1.9549,2.058,2.058];
t=0:1:500;
v=spline(x,y,t); %Geschwindigkeitskurve
%Ableitung der approx. Kurve
a1=diff(v);
a=[a1,0];
t=0:1:500;
Driftwinkel.m
%Driftwinkel für MS Beluga Revolution
%loaded condition, deep water, beta=16°,20°
%
global beta
x=[0,115,244,355,482];
y=(pi/180).*[0,16,20,20,20];
t=0:1:500;
beta=spline(x,y,t);
Geschwindigkeit.m
%Geschwindigkeitverlauf für MS Beluga Revolution
%loaded condition, deep water
%
global v
x=[0,115,244,355,482,500];
y=[6.173,3.7554,3.1381,1.9549,2.058,2.058];
t=0:1:500;
v=spline(x,y,t);
Kruemmung.m
%Krümmungsradiusverlauf für MS Beluga Revolution
%loaded condition, deep water, r approx.236m
%
global r
x=[0,115,244,355,482,500];
y=[2000,300,236,236,236,236];
t=0:1:500;
r=spline(x,y,t);
und schließlich das function-file für den aufrichtenden Hebelarm, in dem das gesuchte [mm] \phi [/mm] ebenfalls mitverwurstet wird: GZ.m
%GZ-Hebelarmkurve
function [h]=GZ(phi)
Schiff;
GM = KM-KG;
BM = KM-0.53*T;
h=(GM + 0.5*BM*(tan(phi)).^2).*sin(phi);
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 23:05 So 20.05.2007 | | Autor: | nschlange |
Wie ist denn die zu erwartende Größenordnung von phi?
Ist das im Bereich bis 130 s realistisch (0,07 Grad)?
Ich würde noch Deine Fktn beta umbenennen, ich hab da beim
ersten Durchlauf immer einen Fehler bekommen, weil es diese
Fktn schon gibt, aber mehr Parameter erwartet.
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 23:16 So 20.05.2007 | | Autor: | nschlange |
Wo kommt überhaupt diese Schleife her?
Ist das dt=1 vielleicht zu groß?
Warum versuchst Du nicht die DGL nicht mit einem implementierten numerischen
Verfahren (Klassiker: Runge-Kutta, 'ode45') zu lösen?
mfg
nschlange
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 19:30 Mo 21.05.2007 | | Autor: | cleophus |
Also die Größenordnung ist mit an Sicherheit grenzender Wahrscheinlichkeit falsch. Unser Dozent hat uns nichts Genaueres verraten, aber aus Erfahrung weiß ich, dass das Schiff sich auf jeden Fall stärker krängt. Ich erwarte eigentlich so 3-15°, wobei es idealerweise weniger als 10° krängen sollte, damit es Stabilitätskriterien erfüllt.
Die Schleife hab ich nach Vorlage selbst zusammengebastelt nachdem uns der Algorithmus vorgegeben wurde, aber könntest du mir kurz erklären, wie ich das ode45-Verfahren anwenden könnte? Würde das trotzdem auf die Unterprogramme zugreifen und wie würde ich das in Matlab schreiben?
Ich werd derweil erstmal deinen Rat befolgen, nschlange, und beta umbenennen. Wobei ich das ja schon merkwürdig finde, dass du da eine Fehlermeldung kriegst... Da muss ich nochmal gucken, ob ich nicht früher schon irgendwas verpeilt habe!
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 23:28 Di 22.05.2007 | | Autor: | nschlange |
Guck mal in die Hilfe, da gibt es einige Beispiele zu den ode*-Befehlen.
Hab allerdings nicht beachtet, dass v und r usw auch Funktionen von t sind.
Ich weiss jetzt gerade nicht ob das geht.
Aber wenn die Schleife vorgegeben war...
Bei mir hatte das Umbenennen von beat nur den Effekt, das beim ersten Durchlauf
keine Fehlermeldung mehr kam. Der Graph blieb der gleiche.
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 19:38 Mo 21.05.2007 | | Autor: | cleophus |
Ich bin auch nochmal deinem Tipp nachgegangen und hab mal dt geändert, aber für alles <1 hab ich die Fehlermeldung ''Subscript indices must either be real positive integers or logicals.'' gekriegt (evtl ein Problem mit (k+1)?) und für alle werte >1 zeichnet er mir eine sehr merkwürdige Kurve, deren Größenordnung leider auch nicht passt... Eigentlich soll da eine Kurve rauskommen, die von 0 auf einen Winkel raufgeht und sich dann um einen mittleren Winkel einpendelt.
Kann mir jemand bestätigen, dass mein Umrechenfaktor [mm] (\pi [/mm] / 180°) richtig ist?
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 23:30 Di 22.05.2007 | | Autor: | nschlange |
dt muss ganzzahlig sein, weil Du ja auf die Komponenten von Vektoren zugreifst.
Die haben nunmal ganzzahlige Indizes.
Vielleicht kann man ja in der Schleife jedes dt in 0.001*dt ändern.
Keine Ahnung, ob das Sinn macht.
Frag doch mal den, der den Algorithmus vorgegeben hat.
|
|
|
|
