¿Descomposición rápida y simple de Helmholtz-Hodge 2D discreta utilizando FFT?

8

Para un protector de pantalla tonto que estoy tratando de desarrollar, me gustaría generar aleatoriamente una matriz 2D de vectores 2D sin divergencia, y luego usarla para generar un gráfico de convolución integral de línea. He escuchado que una forma de hacer esto es generar ruido aleatorio y luego proyectar el componente solenoidal de su descomposición Helmholtz-Hodge. Para hacer eso, intenté usar el siguiente razonamiento:1

Una función tiene descomposición de Helmholtz-Hodge donde y donde son funciones escalares. Por ahora, suponga que el componente armónico desaparece.2 f = h + ϕ + J ψ J = ( 0 - 1 1 0 ) ϕ , ψ hf:R2R22

f=h+ϕ+Jψ
J=(0110)
ϕ,ψh

En el espacio de Fourier, esto se convierte en y podemos definir una proyección solenoidal operador en el espacio de Fourier como que proyecta una función en su componente solenoidal, a través de P = I - kk

Ff=ikϕ^iJkψ^
F-1PFf=Jψ.
P=Ikkkk
F1PFf=Jψ.

Luego intenté implementar esto en Mathematica, aplicándolo a una matriz aleatoria . Primero genero la matriz aleatoria, y aplico la FFT a cada uno de sus dos componentes:21×21×2

arr = RandomReal[{-1, 1}, {2, 21, 21}];
fArr = Fourier /@ arr;

Luego defino en función del índice de matriz:k

k[k1_, k2_] := Mod[{k1 - 1, k2 - 1}, 21, -10]/21;

Luego realizo la proyección en los componentes de Fourier (la singularidad en se deja sola usando una declaración):k=0If

dat = Transpose[
   Table[If[k1 == 1 && k2 == 1, fArr[[;; , k1, k2]], 
     fArr[[;; , k1, k2]] - 
      k[k1, k2] (k[k1, k2].fArr[[;; , k1, k2]])/(k[k1, k2].k[k1, 
           k2])], {k1, 21}, {k2, 21}], {2, 3, 1}];

Entonces iFFT los dos componentes:

projArr = InverseFourier /@ dat;

Esto proporciona una matriz puramente real, y yo ingenuamente esperaría que el resultado fuera una aproximación de . Mi pregunta es:Jψ

  • ¿En qué sentido el resultado se aproxima a ?Jψ

Supuestamente, la descomposición Helmholtz-Hodge de datos 2D no es una tarea trivial, ya que se supone que la rutina HH_DECOMP de Chris Beaumont usa FFT para realizar la descomposición Helmholtz-Hodge, pero también dice (en los comentarios en la parte superior del código) que el método parece incorrecto. Del mismo modo, existen métodos de variación más complicados para realizar descomposiciones de Helmholtz-Hodge de datos 2D, lo que parece sugerir que el método FFT más simple es de alguna manera inadecuado. ¿Por qué? ¿En qué se equivoca el método FFT? ¿Y está mal suponer que el componente armónico desaparece para mi ruido aleatorio?3

(1): Fluidos estables , Jos Stam.

(2): Detección de características en campos vectoriales utilizando la descomposición Helmholtz-Hodge , Alexander Wiebel, página 12.

(3): Descomposición discreta del campo vectorial multiescala , Yiying Tong.

ContenedorDoofus
fuente
Pensaré en tu pregunta real en un minuto, pero ¿por qué no generar una aleatoria y definir ? Lectura adicional: " Ruido de rizo para flujo de fluido de procedimiento ", R. Bridson et al. f = J ψψf=Jψ

Respuestas:

5

La forma más fácil de generar una función 2D libre de divergencia es usar una función de flujo, . Dóndeψ

u=ψyandv=ψx

Una vez que genera un aleatorio , puede usar su representación de Fourier (o diferencias finitas) para calcular las dos derivadas. Sin embargo, no he probado este tipo de cosas con entradas verdaderamente aleatorias. Puede beneficiarse comenzando con una entrada más suave pero aleatoria (como una suma de funciones sinusoidales con frecuencias y coeficientes aleatorios).ψ

Bill Barth
fuente
Gracias, como dices, usando donde es una suma de gaussianos o algo definitivamente funcionaría y sería fácil, y probablemente se vería bien. En cuanto a garantizar la suavidad, una segunda posibilidad es usar donde es un filtro de paso bajo y es ruido aleatorio , que encontré también funciona bien. Una tercera posibilidad que también garantiza condiciones de contorno periódicas es utilizar una suma ponderada de varias de las funciones de la base solenoidal real . JψψF1LPFfLf2Jk^N1N2cos(2πkr+τ)
DumpsterDoofus
4

Hiciste varias preguntas en tu publicación, por lo que si bien Bill ha proporcionado una solución al problema subyacente, creo que alguien también debería decir algunas cosas sobre tus preguntas planteadas.

¿En qué sentido el resultado se aproxima a ?Jψ

Si interpreta la matriz de entrada como el muestreo sin alias de un campo vectorial periódico, ilimitado en banda , entonces el resultado es precisamente el muestreo de . Después de todo, la FFT de la matriz muestreada es esencialmente la serie de Fourier de la función periódica subyacente (suponiendo un muestreo suficiente).h + J ψf=h+ϕ+Jψh+Jψ

Tenga en cuenta que dado que es periódico, es simplemente un campo vectorial constante.hfh

Se supone que la rutina HH_DECOMP de Chris Beaumont usa FFT para realizar la descomposición de Helmholtz-Hodge, pero también dice (en los comentarios en la parte superior del código) que el método parece inexacto.

No estoy seguro de lo que está sucediendo allí porque no puedo leer IDL, pero parece que los campos del vector de prueba no son periódicos.

Del mismo modo, existen métodos variados más complicados para realizar descomposiciones de datos 2D de Helmholtz-Hodge, lo que parece sugerir que el método FFT más simple es de alguna manera inadecuado. ¿Por qué?

Creo que te estás preocupando demasiado. En la segunda página del documento que cita, dicen: "La descomposición discreta de Helmholtz-Hodge en cuadrículas regulares ya se ha utilizado en gráficos (ver [25, 10] por ejemplo) y es relativamente sencillo de implementar con un enfoque de diferencia finita Sin embargo, es mucho más difícil diseñar un método práctico y preciso para cuadrículas arbitrarias ".

¿En qué se equivoca el método FFT?

Si tiene condiciones límite periódicas, nada. De hecho, en esas condiciones, los métodos espectrales como este le brindan una convergencia exponencial en comparación con la convergencia polinómica de los métodos de diferencias finitas. Sin embargo, son más difíciles de aplicar a geometrías arbitrarias.

¿Y está mal suponer que el componente armónico desaparece para mi ruido aleatorio?

Técnicamente, sí, pero en el caso periódico el único campo vectorial armónico posible es una constante. Entonces estás bien dejando solo el componente DC.


fuente