Zufällige 3D Vektoren innerhalb einer Kugel

µ schrieb:

fdfdg schrieb:

krümelkacker schrieb:

Radius ziehen:
Die Wahrscheinlichkeitsdichte für den Radius muss mit dem Radius quadratisch zunehmen, also:
P(r) \propto r^2
da die Kugeloberfläche ja mit dem Radius auch quadratisch zunimmt. Sonst bekommt man hinterher keine Gleichverteilung.

Wie bastelt man das dann im Endeffekt? Wenn Radius zwischen r1 und r2 liegen soll: Kann man einfach gleichverteilt zwischen r1*r2 und r2*r2 ziehen und daraus die Wurzel als Radius nehmen? Müsste doch eine quadratische Verteilung ergeben.

Stell eine Wahrscheinlichkeitsdichte Proportional zu r^2 auf und bestimme den Faktor mit der Bedingung, dass das Integral der Dichte von r1 bis r2 gleich 1 wird. Dann wendet man die Inversionsmethode an.

http://de.wikipedia.org/wiki/Inversionsmethode

Oder kurz: Ja, das stimmt so, fdfdg. Bis auf den Typo und die exakten Exponenten halt

Wen es interessiert, hab hier kurz was gebastelt: http://ideone.com/IsAj2

Werner Salomon

otze schrieb:

Für manche Dinge braucht es nicht C++11 sondern einfache Mathematik.

genau .. und es kommt öfter vor, als man denkt

Folgender Code sollte in etwa das sein, was PI braucht. Er sollte Positionen in einer Kugelschale zwischen zwei gegebenen Radien R1 und R2 erzeugen, die im Volumen annähernd gleich verteilt sind.

Im Gegensatz zu pumuckl und fdfdg bin ich aber der Meinung, dass man beim Radius die 3.Potenz nehmen muss, da der Radius im Volumen mit der 3.Potenz eingeht. Bei einer Kugelschale mit den Innenradius 1 und Außenradius 2 befinden sich 50% des Volumens jenseits eines Radius von ca. 1,651 bzw. = 3.Wurzel( 1 + (2^3 - 1^3)/2 ).

#include <boost/random/uniform_real.hpp>
#include <boost/random/linear_congruential.hpp> // minstd_rand

#include <iostream>
#include <cmath>

class Vector // nur für Demozwecke
{
public:
    typedef double value_type;
    Vector( double x, double y, double z )
        : m_x(x), m_y(y), m_z(z)
    {}
    friend Vector operator*( double factor, const Vector& v )
    {
        return Vector( v.m_x * factor, v.m_y * factor, v.m_z * factor );
    }
    friend double abs( const Vector& v )
    {
        return std::sqrt( v.m_x * v.m_x + v.m_y * v.m_y + v.m_z * v.m_z );
    }
    friend std::ostream& operator<<( std::ostream& out, const Vector& v )
    {
        return out << v.m_x << " " << v.m_y << " " << v.m_z;
    }
private:
    double m_x, m_y, m_z;
};

namespace
{
    const double PI = std::acos(-1.0);
} 

// -- Verteilung zur Generierung eines 'gleichverteilten' Einheitsvektors
template< typename V, typename T = typename V::value_type >
class Sphere
{
public:
    typedef T input_type;
    typedef V result_type;

    Sphere() 
        : m_sin_phi( -1., +1. )
        , m_beta( -PI, PI )
    {}
    template< typename G >
    result_type operator()( G& gen )
    {
        using std::cos;
        using std::sin;
        using std::sqrt;
        // siehe auch http://de.wikipedia.org/wiki/Kugelkalotte
        const T sin_phi = m_sin_phi( gen ); // sin( Breitengrad ) ist gleichverteilt, da H=2*PI*sin(phi)
        const T cos_phi = sqrt( 1 - sin_phi*sin_phi );
        const T beta = m_beta( gen ); // Längenwinkel
        return result_type( cos( beta )*cos_phi, sin( beta )*cos_phi, sin_phi );
    }
private:
    boost::uniform_real< T > m_sin_phi;
    boost::uniform_real< T > m_beta;
};

// --   Verteilung zur Generierung eines Radius im Intervall [x0, x1)
//      mit dem Ziel im Volumen einer Kugelschale eine 'Gleichverteilung' zu erreichen
template< typename T >
class Cubic
{
public:
    typedef T input_type;
    typedef T result_type;
    Cubic( result_type x0, result_type x1  )
        : m_distribution3( x0 * x0 *x0, x1 * x1 * x1 )
    {}
    template< typename G >
    result_type operator()( G& gen )
    {
        using std::pow;
        return pow( m_distribution3( gen ), result_type( 1/3. ) );
    }
private:
    boost::uniform_real< T > m_distribution3;
};

int main()
{
    using namespace std;
    const double R1 = 1.5;
    const double R2 = 4;
    boost::minstd_rand gen;
    Sphere< Vector > sphere;
    Cubic< double > cubic( R1, R2 );
    cout << fixed;
    cout.precision( 2 ); // 2 NK reichen hier
    for( int i=0; i<20; ++i )
    {
        Vector pos = cubic( gen ) * sphere( gen );
        cout << pos << "  " << abs(pos) << endl;
    }
    return 0;
}

Gruß
Werner

@Edit: in den Zeilen 42 und 58 den Typ 'Vector' gegen das Template-Argument V bzw. result_type ersetzt

SeppJ

@Werner: Bist du da sicher, dass die Verteilung deiner Längen stimmt? Dein Code ist mir zu verschlungen, um ihn jetzt noch zu entziffern. Du ziehst nur die dritte Wurzel, oder? Das würde nämlich nicht stimmen, da du dann nämlich entweder die Verteilung bei R1 als 0 hast oder die Verteilung verschoben (und somit die falsche Steigung) hast, je nachdem, was du da genau machst (sorry, wenn ich gerade keine 100 Zeilen Templates dekodieren möchte - habe heute Abend schon genug Denksport betrieben).

Ich habe jetzt stundenlang an der exakten Verteilung der Längen rumgetüftelt, ich kann nicht glauben, dass du so einen einfachen Code heraus bekommst, ich hingegen dieses Monstrum:

#include <iostream>
#include <vector>
#include <functional>
#include <numeric>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <complex>
using namespace std;

int main ()
{
  vector<double> bins(1000,0.);
  const double r1 = 2;
  const double r2 = 5;
  double r1_3 = r1*r1*r1;
  double len = r2-r1;
  double cut = r1*r1 / (pow(r2,3)/3. - r1_3/3.);  // Wert der Verteilung bei r1 -> Offset
   double fac1 = cut * len;   // Der Offset muss noch mit der Länge normiert werden
  for(unsigned i = 0; i < 1000000; ++i)
    {
      double a = rand() *1. / RAND_MAX;
      if (a < fac1)              // Gleichverteilungsanteil der Verteilung
        a = (a/fac1)*len + r1;   // Dies ist noch einfach...
      else                       // ...aber jetzt kommt die x^2 Funktion die da drüber gelegt ist: OMG!
        {
          a = (a-fac1) / (1-fac1); // a ist wieder Einheitsverteilung
          // Herleitung: Mathematica, Zettel und Stift und viel Geduld
          double fac2 = len*len*(2*r1+r2)*a;
          complex<double> fac3 = pow(-2*r1_3+fac2+sqrt(complex<double>(fac2*(-4*r1_3+fac2),0.)),1./3.);
          complex<double> result = 0.5 * (
                                          pow(2,2./3.)*fac3
                                          +2*r1*
                                          (-1. + pow(2,1./3.)*r1 / fac3 )
                                  );
          a = result.real();   // Zusatzwissen: Imaginärteile heben sich genau weg
          a += r1;
        }
      ++bins[a*100];
    }
  double factor = 1./(accumulate(bins.begin(), bins.end(), 0.)*0.01);
  transform (bins.begin(), bins.end(), bins.begin(), bind1st(multiplies<double>(),factor));
  for (unsigned i = 0; i < bins.size(); ++i)
    {
      cout << i / 100. << '\t' << bins[i] << '\n';
    }
}

Dieser Code ist erfolgreich geprüft (beziehungweise er ist mein Prüfcode, einfach mal die Werte plotten, wer es nachvollziehen möchte) gegen die theoretische Wahrscheinlichkeitsverteilung von

P(R) = R^2/(R2^3/3 - R1^3/3) * Theta(R-R1) * Theta(R2-R)

Irgendwie muss man auch die komplexwertige Zwischenrechnung vereinfachen können, aber das würde ich selber gerne sehen und Mathematica schafft es nicht.

Die Längen kann man dann mit Punkten auf einer Einheitssphäre kombinieren, wie hier schon gezeigt wurde. Herleitung findet man hier:
http://mathworld.wolfram.com/SpherePointPicking.html
Das ist zum Glück ein ganzes Ende einfacher.

Ob sich die Rechnung nun lohnt sei dahingestellt. Bei der Wegwerfmethode kann man höchstens (4/3 r^3)/(2r)3 = 1/6 der Zufallszahlen benutzen, wenn man eine Vollkugel hat. Bei einer dünnen Schale sinkt der Anteil sehr schnell nahe 0. Und man hat natürlich auch milde Rechenkosten für die Abstandsberechnung.

Entschuldigung auch an alle Leser, dass mein Code ebenso unleserlich ist wie Werners, wenn auch aus weitaus weniger künstlerisch anspruchsvollen Gründen :p

P.S.: Herleitung ohne funktionierende LaTeX-Tags möchte ich niemandem antun. Im Prinzip ist es nichts kompliziertes, bloß eklig lang ohne CAS (und mit CAS auch ). Die Grundzüge wurden hier im Thread schon erklärt oder Verlinkt, das einzige was glaube ich noch nicht genannt wurde ist, wie man das macht, wenn die Inverse der kumulierten Wahrscheinlichkeitsverteilung bei 0 nicht 0 ist. Das ist das was hier im ersten if-Zweig behandelt wird, indem eine Gleichverteilung passender Höhe "abgeschröpft" und einfach (nach ein bisschen Skalieren) durchgeleitet wird. Der zweite if-Zweig ist dann der Teil der übrig bleibt und mit herkömmlichen Mitteln behandelt werden kann (auch wenn die herkömmlichen Mittel hier ganz schön komplexe Rechnung erfordern). Die Tatsache, dass die Imaginärteile bei negativer Wurzel sich gegenseitig weg heben habe ich zugegebenermaßen nur durch Ausprobieren getestet, aber konstruktionsbedingt muss sich das auch so ergeben, weil die Methodik an sich bewiesen ist und in der Rechnung selbst bis dahin keine Fehler sind (weil der Computer gerechnet hat ).

krümelkacker

SeppJ: Probier's doch mal so zu trenne, wie ich es vorgeschlagen habe. Wenn Du X/Y/Z Gauß-verteilt ziehst und dann normalisierst, hast Du schon auf der Kugeloverfläche eine Gleichverteilung. Dann fehlt nur noch der Abstand zum Ursprung. Ich habe den letzten Schritt nicht zu Ende gerechnet, aber so kompliziert dürfte das nicht sein.

*Rewind*

@Werner: wie lange hast du für den Code gebraucht?

otze

Bashar schrieb:

otze schrieb:

Für manche Dinge braucht e snicht C++11 sondern einfache Mathematik.

"einfache"? Ich hab grad keine Ahnung, wie das mathematisch gehen könnte.

ich meinte in der Tat die ziehen und Wegwerfmethode. Bei den typischen schnellen Zufallszahlengeneratoren ist die analytische Lösung wahrscheinlich langsamer als wegwerfen und neu ziehen. In dem Setup hier ist das auch gut möglich. Ist auch nicht weniger mathematisch, läuft dann halt unter Akzeptanzkriterium.

SeppJ

krümelkacker schrieb:

Ich habe den letzten Schritt nicht zu Ende gerechnet, aber so kompliziert dürfte das nicht sein.

Ja, das hatte ich auch gedacht:

SeppJ schrieb:

hier ist das bei der parabelfoermigen Verteilung der r-Koordinate aber noch einfach

Dankste .

Die Gleichverteilung auf der Kugelfläche ist nämlich der einfache Teil, auch wenn es zuerst nicht so aus sieht.

otze

SeppJ schrieb:

Die Gleichverteilung auf der Kugelfläche ist nämlich der einfache Teil, auch wenn es zuerst nicht so aus sieht.

Wenn wir problemlos für einen gezogenen Radius normalverteilte Punkte auf der Oberflöche erzeugne können, dann müssen wir doch nur den radius ordentlich ziehen.

Mein Ansatz wäre: wir teilen das Volumen der Kugel in unendlich viele Kugeloberflächen auf. dann wird ein Wert des Radius mit dem Verhätlnis gezogen, den seine (infinitesimal dicke) Kugeloberfläche am Integral über alle Oberflächen hat. Das ist praktischerweise das Volumen der Schicht

angenommen der maximale Radius ist 1 und a ist der innere Radius der Kugel, 0 < a < 1

also p(r) = 4* pi * r^2 /(4/3*pi - 4/3*pi *a^3) = 3*r^2/(1-a3), a <= r <= 1

darüber dann hübsch die Inversionmethode anwenden.

Sei q Normalverteilt in U(a,1)

q = \int_a^r' 3r^2/(1-a3) d r
<=> q = (r'^3 - a^3)/(1-a3)
<=> r = (q(1-a^3)+a3)^(1/3)

vermutlich habe ich mich irgendwo verrechnet.

Und dann 2 Winkel gleichverteilt ziehen, zusammen mit dem Radius müsste das einen gleichverteilten Vektor auf der Kugelschicht mit Radius r erzeugen.

Werner Salomon

SeppJ schrieb:

@Werner: Bist du da sicher, dass die Verteilung deiner Längen stimmt? Dein Code ist mir zu verschlungen, um ihn jetzt noch zu entziffern. Du ziehst nur die dritte Wurzel, oder? Das würde nämlich nicht stimmen, da du dann nämlich entweder die Verteilung bei R1 als 0 hast oder die Verteilung verschoben (und somit die falsche Steigung) hast, je nachdem, was du da genau machst

Ja ich ziehe die 3.Wurzel aus einer Verteilung von x0^3 bis x1^3 (s. Cubic<>) - das ist genau das, was otze vorschlägt

otze schrieb:

Mein Ansatz wäre: wir teilen das Volumen der Kugel in unendlich viele Kugeloberflächen auf. dann wird ein Wert des Radius mit dem Verhältnis gezogen, den seine (infinitesimal dicke) Kugeloberfläche am Integral über alle Oberflächen hat. Das ist praktischerweise das Volumen der Schicht

was dann kommt habe ich zwar nicht verstanden - insbesondere die Ausgangsgleichung, aber das Ergebnis

otze schrieb:

<=> r = (q*(1-a^3)+a3)^(1/3)

ist identisch mit dem, was die Verteilung class Cubic<> macht - siehe Zeile 74 und 80. Man ersetze a in otzes Ergebnis durch das Verhältnis x0/x1 aus Cubic<> und rechne es selber nach.

otze schrieb:

vermutlich habe ich mich irgendwo verrechnet.

nö - bis hierher anscheinend nicht.

otze schrieb:

Und dann 2 Winkel gleichverteilt ziehen, zusammen mit dem Radius müsste das einen gleichverteilten Vektor auf der Kugelschicht mit Radius r erzeugen.

warum das falsch ist, ist am Anfang von http://mathworld.wolfram.com/SpherePointPicking.html erklärt und durch eine Graphik demonstriert.

SeppJ schrieb:

Die Längen kann man dann mit Punkten auf einer Einheitssphäre kombinieren, wie hier schon gezeigt wurde. Herleitung findet man hier:
http://mathworld.wolfram.com/SpherePointPicking.html
Das ist zum Glück ein ganzes Ende einfacher.

Danke für den Link SeppJ,
Dort wird genau das beschrieben, was ich in der Klasse Sphere<> codiert habe, nur das die Bezeichner andere sind und die Intervalle der Winkel phasenverschoben sind - z.B. [0,2PI) statt [-PI,+PI). Aber das sollte das Ergebnis nicht verändern.

[Rewind] schrieb:

@Werner: wie lange hast du für den Code gebraucht?

weiß ich nicht mehr genau. Weniger als eine Stunde auf jeden Fall - vielleicht 45min, die meiste Zeit habe ich dafür gebraucht, aus der Formel für die Oberfläche der Kugelkalotte die Verteilung für den Einheitsvektor zu zaubern. Das ist das, was in den Zeilen 55 bis 57 steht.

Gruß
Werner

Werner Salomon

SeppJ schrieb:

Ich habe jetzt stundenlang an der exakten Verteilung der Längen rumgetüftelt, ich kann nicht glauben, dass du so einen einfachen Code heraus bekommst, ich hingegen dieses Monstrum:

Ich habe noch mal ein kleines Testprogramm geschrieben.
Es ist wirklich so, dass Dein 'Monstrum' exakt die selbe Verteilung liefert wie meine Verteilung class Cubic<>. Bei der identischen Zufallszahl liefern unsere Algorithmen zwar unterschiedliche Radien, aber bei näherer Betrachtung muss das auch so sein, aber die Verteilung ist die selbe.
Also Du hast Dich nicht verrechnet.

SeppJ schrieb:

Die Grundzüge wurden hier im Thread schon erklärt oder Verlinkt, ...

Das kann ich in keiner Weise nachvollziehen. Kannst Du da noch ein wenig konkreter werden, wie die Ausgangsformel oder Idee ist, die dann zu Deiner Lösung führt. Das würde mich interessieren.

Gruß
Werner

SeppJ

Werner Salomon schrieb:

SeppJ schrieb:

Die Grundzüge wurden hier im Thread schon erklärt oder Verlinkt, ...

Das kann ich in keiner Weise nachvollziehen. Kannst Du da noch ein wenig konkreter werden, wie die Ausgangsformel oder Idee ist, die dann zu Deiner Lösung führt. Das würde mich interessieren.

Im Nachhinein gesehen war's eher umständlich, auf otzes Trick hätte ich kommen müssen. Mein's ist im Prinzip die Inversionsmethode direkt mit dem Kopf voran durch die Betonwand gerechnet, anstatt die Tür zu nehmen .

Die Verteilung die ich haben möchte ist

P(R) = R^2/(R2^3/3 - R1^3/3) * Theta(R-R1) * Theta(R2-R)

Ich habe eine GLeichverteilung 0-1 die ich dazu umformen möchte. Wenn ich nun für die Inversionsmethode die kumulierte Wahrscheinlichkeitsverteilung zu berechnen, bekomme ich eine Funktion die bis R1 0 ist und dann wie x^3 bis 1 bei R2 steigt. Leider nicht direkt umkehrbar, weil mehrdeutig bei 0.

Und an dieser Stelle hätte ich es wie otze machen sollen . Was ich stattdessen gemacht habe: Das Problem, wie ich es sah, war, dass die Verteilung bei R1 nicht 0 ist. Wenn's einfach nur quadratisch von R1 bis R2 ginge und P(R1)=0, wäre das Problem einfach. Aber es ist im Prinzip ein "Block" von Gleichverteilung unter der quadratischen Verteilung. Daher habe ich gesagt, ich nehme aus meinen gleichverteilten Zufallszahlen einen Block eben dieser Größe (nämlich P(R1)*(R2-R1)) heraus und strecke diesen, so dass er diesem Block unter meiner Parabel entspricht. Das ist der erste if-Zweig bei mir.

Den Rest der gleichverteilten Zufallszahlen habe ich dann wieder gestreckt, so dass ich wieder eine Gleichvertelung habe, so dass ich die Inversionsmethode nach Schulbuch anwenden kann (Zeile 27). Nun hatte ich nur noch den parabelförmigen Teil von P(R) um den ich mich kümmern musste. Auf den Nullpjunkt verschoben und die Inversionsmethode direkt angewendet. Doch dummerweise haben sich in der Rechnung bis dahin derart viele Vorfaktoren angesammelt, dass das dann doch nicht mehr so einfach war, wie der Plan vorsah. Als kumulierte Wahrscheinlichkeitsfunktion erhielt ich

F(R) = (r^2 (r + 3 R1))/((R1 - R2)^2 (2 R1 + R2))

Nun, und das habe ich dann umgekehrt . Und während draußen Gewitterstürme tobten und der strenge Geruch von Formaldehyd von den eingemachten Organen im Regal aufstieg, erwachte so in meinem Labor beim kalten Lichte der Gasüberschlagslampen das MONSTER . Muahuahuahuahua!!! Sie haben mich wahnsinnig genannt!! Ich werde ihnen zeigen, wer hier wahnsinnig ist!!!!!

Werner Salomon

SeppJ schrieb:

Die Verteilung die ich haben möchte ist
P(R) = R^2/(R2^3/3 - R1^3/3) * Theta(R-R1) * Theta(R2-R)

wo kommt diese Gleichung her. Und was ist die Funktion Theta()? In Statistik bin ich leider ziemlich schwach.

Gruß
Werner

SeppJ

Heaviside-Stufenfunktion .

Die Verteilungsfunktion ergibt sich dadurch, dass die Kugelflächen wie R^2 gehen, die Verteilung der Radien die wir ziehen, muss daher auch proportional zu R^2 sein, damit man hinterher eine gleichmäßige Verteilung im Raum bekommt. Die Werte sollen durch R1 und R2 eingeschränkt sein, daher die Stufenfunktionen. Und die Gesamtwahrscheinlichkeit muss 1 sein, das ist die Normierung unter dem Bruch.

*Rewind*

SeppJ schrieb:

Nun, und das habe ich dann umgekehrt . Und während draußen Gewitterstürme tobten und der strenge Geruch von Formaldehyd von den eingemachten Organen im Regal aufstieg, erwachte so in meinem Labor beim kalten Lichte der Gasüberschlagslampen das MONSTER . Muahuahuahuahua!!! Sie haben mich wahnsinnig genannt!! Ich werde ihnen zeigen, wer hier wahnsinnig ist!!!!!

Zu viel Sinn City geguckt?

volkard

SeppJ schrieb:

Ob sich die Rechnung nun lohnt sei dahingestellt. Bei der Wegwerfmethode kann man höchstens (4/3 r^3)/(2r)3 = 1/6 der Zufallszahlen benutzen, wenn man eine Vollkugel hat.

Da ist Π verlorengegangen.
(4/3 Π r^3)/(2r)3 = 1/6 Π = 52%

Im Fall, daß die innere Kugel den halben Radius hat, sind's noch 46%, wie pumuckl auf Seite 2 festgestellt hatte.

otze

Werner Salomon schrieb:

was dann kommt habe ich zwar nicht verstanden - insbesondere die Ausgangsgleichung

Okay, ich geb mir mal mehr Muehe

Eine vereinfachende Definition zuerst, ich bezeichne die Oberflaeche einer Kugel mit Radius r mit:
A(r) = 4* pi * r^2

was ich nun machen will ist, wie du bereits sagtest: die Auftrittswahrscheinlichkeit eines r muss proportional zu A(r) sein. Da wir zugleich die Nebenbedingung haben, dass das \int_a^1 p(r) dr = 1 sein muss, ergibt sich der Proportionalitaetsfaktor als die Normalisierungskonstante.

p(r) = A(r) / \int_a^1 A(r) dr = A(r)/ (4/3*pi - 4/3*pi*a^3)

und dann nurnoch A(r) einsetzen.

Jetzt bin ich aber auf deinen Gedankengang gespannt. Du scheinst irgendwie straight forward gewusst zu haben, wie die r Verteilt sein muessen.

brotbernd

volkard schrieb:

SeppJ schrieb:

Ob sich die Rechnung nun lohnt sei dahingestellt. Bei der Wegwerfmethode kann man höchstens (4/3 r^3)/(2r)3 = 1/6 der Zufallszahlen benutzen, wenn man eine Vollkugel hat.

Da ist Π verlorengegangen.
(4/3 Π r^3)/(2r)3 = 1/6 Π = 52%

Im Fall, daß die innere Kugel den halben Radius hat, sind's noch 46%, wie pumuckl auf Seite 2 festgestellt hatte.

Zieht man in einer Würfelschale, die die Kugelschale enthält sind es auch wieder 52%. Man kann also auch bei sehr dünnen Schalen genauso viele Zufallszahlen nutzen.

Dobi

Wenn die Kugelschale sehr dünn wird, wird man auch bei der Würfelschalenmethode ziemlich viele Zahlen wegwerfen müssen. Man könnte so eine dünne Kugelschale dann mit Würfeln approximieren. Ist dann dann aber vermutlich kaum noch einfacher als die hier diskutierte elegante/komplizierte Lösung.

brotbernd

Dobi schrieb:

Wenn die Kugelschale sehr dünn wird, wird man auch bei der Würfelschalenmethode ziemlich viele Zahlen wegwerfen müssen. Man könnte so eine dünne Kugelschale dann mit Würfeln approximieren. Ist dann dann aber vermutlich kaum noch einfacher als die hier diskutierte elegante/komplizierte Lösung.

Es ist immer 0,52.
4/3 * pi * (r1^3 - r2^3) / (8 * (r1^3 - r2^3) = 1/6 pi

Edit: Falscher Fehler. Der innere Würfel ist natürlich innerhalb der inneren Kugel. Das war also quatsch...

brotbernd

...