1

ある距離 b にある球体と無限円柱の交差部分の体積を計算したかったのですが、手早く汚い python スクリプトを使用してそれを行うことにしました。私の要件は、有効数字が 3 桁を超える 1 秒未満の計算です。

私の考えは次のとおりです。半径 R の球を、その中心が原点にあるように配置し、半径 R' の円柱を、その軸が (b,0,0 ）。円柱の内側にある場合は 1 を返し、そうでない場合は 0 を返すステップ関数を使用して、球を積分します。したがって、球と円柱の両方の内側にあることによって制約されたセット、つまり交差で 1 を積分します。

scipy.intigrate.tplquad を使用してこれを試しました。うまくいきませんでした。次のような警告が表示されるため、ステップ関数の不連続が原因だと思います。もちろん、私はこれを間違っているかもしれません。ばかげた間違いを犯していないと仮定すると、交差点の範囲を定式化して、ステップ関数の必要性を取り除くことができますが、最初にフィードバックを得ようとするかもしれないと考えました. 誰でも間違いを見つけたり、簡単な解決策を指摘したりできますか?

警告: サブディビジョンの最大数 (50) に達しました。
制限を増やしても改善が得られない場合は
、問題を判断するために被積分関数を分析することをお勧めします。局所的な困難の位置 (特異点、不連続点) を決定できる場合は、区間を分割し、部分範囲で積分器を呼び出すことでおそらく利益が得られます。おそらく、専用のインテグレーターを使用する必要があります。

コード：

from scipy.integrate import tplquad
from math import sqrt


def integrand(z, y, x):
    if Rprim >= (x - b)**2 + y**2:
        return 1.
    else:
        return 0.

def integral():
    return tplquad(integrand, -R, R, 
                   lambda x: -sqrt(R**2 - x**2),          # lower y
                   lambda x: sqrt(R**2 - x**2),           # upper y
                   lambda x,y: -sqrt(R**2 - x**2 - y**2), # lower z
                   lambda x,y: sqrt(R**2 - x**2 - y**2), # upper z
                   epsabs=1.e-01, epsrel=1.e-01
                   )

R=1
Rprim=1
b=0.5
print integral()
4

4 に答える 4

3

[0, 0, 0]球の原点が にあり、その半径がであるような方法でデータを変換およびスケーリングできると仮定すると1、単純な確率的近似で合理的な答えが十分に速く得られる可能性があります。したがって、線に沿った何かが良い出発点になる可能性があります。

import numpy as np

def in_sphere(p, r= 1.):
    return np.sqrt((p** 2).sum(0))<= r

def in_cylinder(p, c, r= 1.):
    m= np.mean(c, 1)[:, None]
    pm= p- m
    d= np.diff(c- m)
    d= d/ np.sqrt(d** 2).sum()
    pp= np.dot(np.dot(d, d.T), pm)
    return np.sqrt(((pp- pm)** 2).sum(0))<= r

def in_sac(p, c, r_c):
    return np.logical_and(in_sphere(p), in_cylinder(p, c, r_c))

if __name__ == '__main__':
    n, c= 1e6, [[0, 1], [0, 1], [0, 1]]
    p= 2* np.random.rand(3, n)- 2
    print (in_sac(p, c, 1).sum()/ n)* 2** 3
于 2011-04-26T15:36:55.933 に答える
2

2 つのドメインにわたって一定である不連続関数に対してトリプル アダプティブ数値積分を実行することは、特に速度または精度のいずれかを確認したい場合は、非常に悪い考えです。

問題を分析的に削減することは、はるかに優れたアイデアだと思います。

変換により、円柱を軸に合わせます。これにより、球が原点以外の点に移動します。

次に、その軸に沿って球と円柱が交差する限界を見つけます。

その軸変数を積分します。軸に沿った任意の固定値での交点の面積は、単純に 2 つの円の交点の面積であり、三角法と少しの努力を使用して簡単に計算できます。

最終的には、計算時間をほとんど必要とせずに、正確な結果が得られます。

于 2011-04-26T14:39:48.357 に答える
1

eat で提案されているように、単純な MC 統合を使用して解決しましたが、実装が遅くなりました。私の要件は増加しました。したがって、ウッドチップが示唆するように、問題を数学的に再定式化しました。

基本的に、x の限界を z と y の関数として定式化し、y を z の関数として定式化しました。次に、本質的に、境界を使用して、交差点で f(z,y,z)=1 を統合しました。これを行ったのは、速度が向上したためで、ボリュームと b をプロットできるようになり、比較的小さな変更でより複雑な機能を統合できるようになったからです。

誰かが興味を持っている場合に備えて、コードを含めます。

from scipy.integrate import quad
from math import sqrt
from math import pi

def x_max(y,r):
    return sqrt(r**2-y**2)

def x_min(y,r):
    return max(-sqrt(r**2 - y**2), -sqrt(R**2 - y**2) + b) 

def y_max(r):
    if (R<b and b-R<r) or (R>b and b-R>r):
        return sqrt( R**2 - (R**2-r**2+b**2)**2/(4.*b**2) )
    elif r+R<b:
        return 0.
    else: #r+b<R
        return r

def z_max():
    if R>b:
        return R
    else:
        return sqrt(2.*b*R - b**2) 

def delta_x(y, r):
    return  x_max(y,r) - x_min(y,r)

def int_xy(z):
    r = sqrt(R**2 - z**2)
    return quad(delta_x, 0., y_max(r), args=(r))

def int_xyz():
    return quad(lambda z: int_xy(z)[0], 0., z_max())

R=1.
Rprim=1.
b=0.5
print 4*int_xyz()[0]
于 2011-04-27T16:16:10.217 に答える
0

まず、交差点の体積を手で計算できます。それをしたくない（またはできない）場合は、次の方法があります。

ドメインの四面体メッシュを生成し、セル ボリュームを合計します。pygalmeshとmeshplexの例(どちらも私が作成):

import pygalmesh
import meshplex
import numpy

ball = pygalmesh.Ball([0, 0, 0], 1.0)
cyl = pygalmesh.Cylinder(-1, 1, 0.7, 0.1)
u = pygalmesh.Intersection([ball, cyl])

mesh = pygalmesh.generate_mesh(u, cell_size=0.05, edge_size=0.1)

points = mesh.points
cells = mesh.cells["tetra"]

# kick out unused vertices
uvertices, uidx = numpy.unique(cells, return_inverse=True)
cells = uidx.reshape(cells.shape)
points = points[uvertices]

mp = meshplex.MeshTetra(points, cells)
print(sum(mp.cell_volumes))

これにより、

ここに画像の説明を入力

2.6567890958740463ボリュームとして印刷されます。精度を上げるには、セルまたはエッジのサイズを小さくします。

于 2017-03-24T13:02:32.063 に答える