Lösung einer DGL stark vom Startwert abhängig

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Axel

Forum-Anfänger


	Beiträge: 12

	Anmeldedatum: 29.07.11

	Wohnort: ---

	Version: ---

Verfasst am: 25.02.2013, 14:53 Titel: Lösung einer DGL stark vom Startwert abhängig

Hallo,

ich möchte mittels folgender DGL die Form einer Blase/eines Tropfen lösen. Die Gleichung dazu ist in der Literatur häugifer zu finden. z.B. in (Meyers, T. G.; Charpin, J. P. F. (2009): A mathematical model of the Leidenfrost effect on an axisymmetric droplet. In: Phys. Fluids 21 (6), S. 063101063101-8.)

Die Gleichung lautet:

$\frac{\partial}{\partial\xi}\begin{bmatrix} \Theta \\ r \\ z \end{bmatrix} = \begin{bmatrix} -\frac{1}{a^2} z + z_0 - \frac{tan(\Theta)}{r sqrt{1+tan^2(\Theta)}}\\ cos(\Theta) \\ sin(\Theta) \end{bmatrix}$

a ist eine Konstante (Stoffwert) und $z_0$ ein zu setzender Vorgabewert, welcher (nach der Theorie) die Größe der Blase bestimmt. $\xi$ beschreibt die Bogenlänge und $\Theta$ den aktuellen Winkel zur r Achse. Es liegt ein r,z Koordinatensystem vor.

Ich versuche diese Gleichung wie folgt zu lösen:

Hauptprogramm: (zur Anschauung etwas reduziert)

Code:

% Bestimmung der Blasenform nach Myers und Charpin
clear all
clc
close all

%% Eingabewerte
rho_l = 1419.82; %Dichte Fluessigkeit
g = 9.81; %Schwerkraft
sigma = 0.008273;

global theta;
y0 = [0; 1e-15; 0];
theta1 = 50; %Randwinkel

clear y
theta = 180-theta1;

%% Startwerte
zetaSpan=[0,1]; %Lösungsbereich Kurvenparameter
a=sqrt(sigma/(rho_l*g)); %Berechnung Kapillarlänge

z2 = 0;
for z_0 = -3333:-1:-3340

%% Berechnung der Tropfengeometrie mit vorgebenem Volumen
z2 = z2+1;

ode=@(zeta,y) Tropfenform(zeta,y,a,z_0);
options=odeset('Events',@events,'Refine',100);%,'InitialStep',1e-6,'MaxStep',1e-4);
[zeta,y]=ode45(ode,zetaSpan,y0,options);

figure1 = subplot(3,2,z2);
hold on
plot(y(:,2)*1e3,(y(:,3))*1e3,'LineWidth',2);
axis equal;
xlim([0 1.5]);
ylim([-2 0]);
text(1,-1,['z_0 = ',num2str(z_0)]);
drawnow

end

Funktion ohne Link?

Tropfenform:

Code:

function dydzeta = Tropfenform(zeta, y, a, z_0)

dydzeta=[-(1/(a^2))*(y(3)) + z_0 - (1/y(2))*tan(y(1))/(sqrt(1+(tan(y(1)))^2)); cos(y(1)); sin(y(1))];

Funktion ohne Link?

events:

Code:

function [value,isterminal,direction] = events(zeta,y)

global theta;

%value=y(1)+pi;
value = y(1)+theta*pi/180;
isterminal =1;
direction=0;

Funktion ohne Link?

In events wird die Abbruchbedingung gesetzt. Ich habe gerade leider kein Webspace um ein bild hochzuladen. Deshalb habe ich das subfigure eingefügt. Mein Problem ist nun, dass ab einem bestimmten wert sich die Lösung strak ändert. Für z_0 <= -3336 erhalte ich die Lösung in der gewünschten Form (eine halbe Blase). Für z_0 > -3336 geht die Funktion für wachsende $\xi$ dann aber leider "in die andere Richtung"...
Ich habe leider keine Ahnung warum, da für diese Werte von z_0 die Blase eigentlich einfach größer werden sollte.

Hat jemand eine Idee, woran dieses liegt oder wie man es umgehen kann?

Schonmal vielen Dank
Gruß

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 25.02.2013, 15:38 Titel:

Hallo,

die DGL wird mit einer gewissen Genauigkeit simuliert. Bei manchen z0 führt der Genauigkeitsverlust dazu, dass du in eine falsche Lösung läufst (kann man sich vermutlich anhand eines Richtungsfelds schön ansehen, z.B. mit quiver).

Abhilfe:

Code:

options=odeset('Events',@events,'Refine',100, 'RelTol', 1e-12 );

Funktion ohne Link?

Grüße,
Harald

Axel

Themenstarter

Forum-Anfänger


	Beiträge: 12

	Anmeldedatum: 29.07.11

	Wohnort: ---

	Version: ---

Verfasst am: 25.02.2013, 16:34 Titel:

Hallo Harald,

danke für die schnelle Antwort. Es verbessert das Problem wirklich. Leider tritt der Fehler aber immernoch auf. Jetzt zwischen -3332 und -3333. Beim weiteren Verkleinern der RelTol bekomme ich dann irgendwann eine Warnung, dass diese auf 2.22045e-014 erhöht wurde.

Danke und Gruß
Axel

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 25.02.2013, 16:57 Titel:

Hallo,

wenn du 'abstol' auch noch verkleinerst, bekommst du auch noch für -3332 die gewünschte Lösung, allerdings für -3330 nicht mehr. Ich vermute, dass man da einfach an die Grenzen der Solver stößt.

Grüße,
Harald

Axel

Themenstarter

Forum-Anfänger


	Beiträge: 12

	Anmeldedatum: 29.07.11

	Wohnort: ---

	Version: ---

Verfasst am: 25.02.2013, 19:02 Titel:

Fällt jemandem evtl. ein anderer Ansatz ein um die Gleichungen zu lösen?

Danke und Gruß
Axel

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 25.02.2013, 22:31 Titel:

Hallo,

was mir noch aufgefallen ist: du verwendest in der Anfangsbedingung 1e-15. Hat das einen tieferen Hintergrund oder ist das einfach "ein Wert nah bei 0"? Falls letzteres, könnte man auch versuchen, diesen Wert zu variieren - etwas näher nach 0 hin oder von 0 weg.

Grüße,
Harald

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 26.02.2013, 11:36 Titel: Re: Lösung einer DGL stark vom Startwert abhängig

Hallo Axel,

Nur eine Randbemerkung: Der leider übliche brutale Clear-Header "clear all, clc, close all" ist kein Fehler, aber hilfreich ist er im Allgemeinen auch nicht. Das Löschen aller Funktionen aus dem Speicher raubt nur Zeit, da das Nachladen von der Platte langsam ist.
Einige Zeilen später dann nochmal "clear y" anzuwenden, bringt dann gar nichts, da ja bereits alles gecleart wurde.

Gruß, Jan

Axel

Themenstarter

Forum-Anfänger


	Beiträge: 12

	Anmeldedatum: 29.07.11

	Wohnort: ---

	Version: ---

Verfasst am: 27.02.2013, 10:08 Titel:

Harald hat Folgendes geschrieben:

Hallo,

du verwendest in der Anfangsbedingung 1e-15. Hat das einen tieferen Hintergrund oder ist das einfach "ein Wert nah bei 0"? Falls letzteres, könnte man auch versuchen, diesen Wert zu variieren - etwas näher nach 0 hin oder von 0 weg.

Das ist einfach ein Wert nahe 0 da [0; 0; 0] eine triviale Lösung bietet (eine Ebene). Ein erhöhen des Startwertes führt irgendwann zu einer Lösung. Dazu muss ich den Wert aber so weit erhöhen, dass die Lösung nicht mehr dem eigentlichen Problem entspticht.

Jan S hat Folgendes geschrieben:

Der leider übliche brutale Clear-Header "clear all, clc, close all" ist kein Fehler, aber hilfreich ist er im Allgemeinen auch nicht. Das Löschen aller Funktionen aus dem Speicher raubt nur Zeit, da das Nachladen von der Platte langsam ist.
Einige Zeilen später dann nochmal "clear y" anzuwenden, bringt dann gar nichts, da ja bereits alles gecleart wurde.

Das mit dem Header war mir nicht bewusst. Danke für den Tipp. Das

Code:

clear y

Funktion ohne Link?

ist noch aus dem kürzen auf das Beispiel übrig geblieben, da ich vorhatte, das Programm mittels for Schleife über mehrere Kontaktwinkel auszuführen. Über

Code:

size(y,1)

Funktion ohne Link?

konnte ich dann gut prüfen, welche Lösung vorlag.


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 - 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.