Resolver la ecuación de Laplace

13

Introducción a la matemática numérica

Este es el "¡Hola Mundo!" de PDE (ecuaciones diferenciales parciales). La ecuación de Laplace o difusión aparece a menudo en física, por ejemplo, ecuación de calor, deformación, dinámica de fluidos, etc. Como la vida real es 3D, pero queremos decir "¡Hola, mundo!" y no cantar "99 botellas de cerveza, ..." esta tarea se da en 1D. Puede interpretar esto como una túnica de goma atada a una pared en ambos extremos con cierta fuerza aplicada.

En un [0,1]dominio, encuentre una función upara la función fuente dada fy los valores límite u_Ly u_Rtal que:

  • -u'' = f
  • u(0) = u_L
  • u(1) = u_R

u'' denota la segunda derivada de u

Esto puede resolverse puramente teórico, pero su tarea es resolverlo numéricamente en un dominio discretizado x para Npuntos:

  • x = {i/(N-1) | i=0..N-1}o basado en 1:{(i-1)/(N-1) | i=1..N}
  • h = 1/(N-1) es el espacio

Entrada

  • f como función o expresión o cadena
  • u_L, u_Rcomo valores de coma flotante
  • N como entero> = 2

Salida

  • Matriz, Lista, algún tipo de cadena separada de utal manera queu_i == u(x_i)

Ejemplos

Ejemplo 1

Entrada: f = -2, u_L = u_R = 0, N = 10(no tome f=-2mal, no es un valor, sino una función constante que los retornos -2para todos xEs como una fuerza de la gravedad constante en nuestra cuerda.).

Salida: [-0.0, -0.09876543209876543, -0.1728395061728395, -0.22222222222222224, -0.24691358024691357, -0.24691358024691357, -0.22222222222222224, -0.1728395061728395, -0.09876543209876547, -0.0]

Existe una solución exacta fácil: u = -x*(1-x)

Ejemplo 2

Entrada: f = 10*x, u_L = 0 u_R = 1, N = 15(Aquí hay una gran cantidad de viento arriba a la derecha)

Salida: [ 0., 0.1898688, 0.37609329, 0.55502915, 0.72303207, 0.87645773, 1.01166181, 1.125, 1.21282799, 1.27150146, 1.29737609, 1.28680758, 1.2361516, 1.14176385, 1.]

La solución exacta para esto dice: u = 1/3*(8*x-5*x^3)

Ejemplo 3

Entrada: f = sin(2*pi*x), u_L = u_R = 1, N = 20(Alguien rompió la gravedad o hay una especie de arriba y la dirección del viento)

Salida: [ 1., 1.0083001, 1.01570075, 1.02139999, 1.0247802, 1.0254751, 1.02340937, 1.01880687, 1.01216636, 1.00420743, 0.99579257, 0.98783364, 0.98119313, 0.97659063, 0.9745249, 0.9752198, 0.97860001, 0.98429925, 0.9916999, 1.]

Aquí la solución exacta es u = (sin(2*π*x))/(4*π^2)+1

Ejemplo 4

Entrada: f = exp(x^2), u_L = u_R = 0,N=30

Salida: [ 0. 0.02021032 0.03923016 0.05705528 0.07367854 0.0890899 0.10327633 0.11622169 0.12790665 0.13830853 0.14740113 0.15515453 0.16153488 0.1665041 0.17001962 0.172034 0.17249459 0.17134303 0.16851482 0.1639387 0.15753606 0.1492202 0.13889553 0.12645668 0.11178744 0.09475961 0.07523169 0.05304738 0.02803389 0. ]

Tenga en cuenta la ligera asimetría

FDM

Un posible método para resolver esto es el método de diferencia finita :

  • reescribir -u_i'' = f_icomo
  • (-u_{i-1} + 2u_i - u{i+1})/h² = f_i que es igual
  • -u_{i-1} + 2u_i - u{i+1} = h²f_i
  • Configura las ecuaciones:

  • Que son iguales a una ecuación matriz-vector:

  • Resuelve esta ecuación y genera el u_i

Una implementación de esto para demostración en Python:

import matplotlib.pyplot as plt
import numpy as np
def laplace(f, uL, uR, N):
 h = 1./(N-1)
 x = [i*h for i in range(N)]

 A = np.zeros((N,N))
 b = np.zeros((N,))

 A[0,0] = 1
 b[0] = uL

 for i in range(1,N-1):
  A[i,i-1] = -1
  A[i,i]   =  2
  A[i,i+1] = -1
  b[i]     = h**2*f(x[i])

 A[N-1,N-1] = 1
 b[N-1]     = uR

 u = np.linalg.solve(A,b)

 plt.plot(x,u,'*-')
 plt.show()

 return u

print laplace(lambda x:-2, 0, 0, 10)
print laplace(lambda x:10*x, 0, 1, 15)
print laplace(lambda x:np.sin(2*np.pi*x), 1, 1, 20)

Implementación alternativa sin Matrix Algebra (usando el método Jacobi )

def laplace(f, uL, uR, N):
 h=1./(N-1)
 b=[f(i*h)*h*h for i in range(N)]
 b[0],b[-1]=uL,uR
 u = [0]*N

 def residual():
  return np.sqrt(sum(r*r for r in[b[i] + u[i-1] - 2*u[i] + u[i+1] for i in range(1,N-1)]))

 def jacobi():
  return [uL] + [0.5*(b[i] + u[i-1] + u[i+1]) for i in range(1,N-1)] + [uR]

 while residual() > 1e-6:
  u = jacobi()

 return u

Sin embargo, puede usar cualquier otro método para resolver la ecuación de Laplace. Si usa un método iterativo, debe iterar hasta el residual |b-Au|<1e-6, bsiendo el vector del lado derechou_L,f_1h²,f_2h²,...

Notas

Dependiendo de su método de solución, es posible que no resuelva los ejemplos exactamente para las soluciones dadas. Al menos para N->infinityel error debe acercarse a cero.

Las lagunas estándar no están permitidas , se permiten incorporaciones para PDE.

Prima

Una bonificación de -30% para mostrar la solución, ya sea gráfica o ASCII-art.

Victorioso

¡Esto es codegolf, por lo que gana el código más corto en bytes!

Karl Napf
fuente
Recomiendo agregar un ejemplo que no sea analíticamente solucionable, por ejemplo, con f(x) = exp(x^2).
falla
@flawr Claro, tiene una solución, sin embargo, la función de error está involucrada.
Karl Napf
1
Lo sentimos, esa fue quizás la expresión incorrecta, ¿podría ser mejor "antiderivada no elemental"? Me refiero a funciones como log(log(x))o sqrt(1-x^4)que tienen una integral, que sin embargo no se puede expresar en funciones elementales.
falla
@flawr No, está bien, la función de error no es elemental, solo quería decir que hay una manera de expresar la solución analíticamente pero u(x) = 1/2 (-sqrt(π) x erfi(x)+sqrt(π) erfi(1) x+e^(x^2)-e x+x-1)no es exactamente calculable.
Karl Napf
¿Por qué iterar hasta 1e-6 y no iterar hasta 1e-30?
RosLuP

Respuestas:

4

Mathematica, 52.5 bytes (= 75 * (1 - 30%))

+0.7 bytes por comentario de @flawr.

ListPlot[{#,u@#}&/@Subdivide@#4/.NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]]&

Esto traza la salida.

p.ej

ListPlot[ ... ]&[10 x, 0, 1, 15]

ingrese la descripción de la imagen aquí

Explicación

NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]

Resolver para la función u.

Subdivide@#4

Subdivide el intervalo [0,1] en N (cuarta entrada) partes.

{#,u@#}&/@ ...

Mapa ua la salida de Subdivide.

ListPlot[ ... ]

Trazar el resultado final.

Solución no gráfica: 58 bytes

u/@Subdivide@#4/.NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]&
JungHwan Min
fuente
Esto no funciona paraf(x) = exp(x^2)
flawr
Quizás desee utilizarlo NDSolvepara el caso general de soluciones no elementales.
error
6

Matlab, 84, 81.2 79.1 bytes = 113 - 30%

function u=l(f,N,a,b);A=toeplitz([2,-1,(3:N)*0]);A([1,2,end-[1,0]])=eye(2);u=[a,f((1:N-2)/N)*(N-1)^2,b]/A;plot(u)

Tenga en cuenta que en este ejemplo utilizo vectores de fila, esto significa que la matriz Ase transpone. fse toma como un manejador de funciones, a,bson las contracciones de Dirichlet del lado izquierdo / derecho.

function u=l(f,N,a,b);
A=toeplitz([2,-1,(3:N)*0]);       % use the "toeplitz" builtin to generate the matrix
A([1,2,end-[1,0]])=eye(2);        % adjust first and last column of matrix
u=[a,f((1:N-2)/N)*(N-1)^2,b]/A;   % build right hand side (as row vector) and right mu
plot(u)                           % plot the solution

Para el ejemplo, f = 10*x, u_L = 0 u_R = 1, N = 15esto resulta en:

falla
fuente
3

R, 123,2 102,9 98,7 bytes (141-30%)

Editar: ¡Ahorré un puñado de bytes gracias a @Angs!

Si alguien quiere editar las imágenes, siéntase libre de hacerlo. Esto es básicamente una adaptación R de las versiones de Matlab y Python publicadas.

function(f,a,b,N){n=N-1;x=1:N/n;A=toeplitz(c(2,-1,rep(0,N-2)));A[1,1:2]=1:0;A[N,n:N]=0:1;y=n^-2*sapply(x,f);y[1]=a;y[N]=b;plot(x,solve(A,y))}

Ungolfed y explicado:

u=function(f,a,b,N){
    n=N-1;                                              # Alias for N-1
    x=0:n/n;                                            # Generate the x-axis
    A=toeplitz(c(2,-1,rep(0,N-2)));                     # Generate the A-matrix
    A[1,1:2]=1:0;                                       # Replace first row--
    A[N,n:N]=0:1;                                       # Replace last row
    y=n^-2*sapply(x,f)                                  # Generate h^2*f(x)
    y[1]=a;y[N]=b;                                      # Replace first and last elements with uL(a) and uR(b)
    plot(x,solve(A,y))                                  # Solve the matrix system A*x=y for x and plot against x 
}

Ejemplo y casos de prueba:

La función nombrada y no reflejada se puede llamar usando:

u(function(x)-2,0,0,10)
u(function(x)x*10,0,1,15)
u(function(x)sin(2*pi*x),1,1,20)
u(function(x)x^2,0,0,30)

Tenga en cuenta que el fargumento es una función R.

Para ejecutar la versión de golf simplemente use (function(...){...})(args)

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

Billywob
fuente
Creo que puede dejar el is.numeric(f)cheque si declara fcomo función, no es necesario pasarlo directamente en la llamada de función al solucionador.
Karl Napf
Ah, ya veo, no sabía que hay una diferencia entre esos dos. Bueno, si es más corto, puede modificar su solucionador para aceptarlo fcomo una función para que no tenga que verificar el caso, es una constante (función).
Karl Napf
1
@Billywob no hay necesidad de fser numérico. f = (function(x)-2)funciona para el primer ejemplo, por lo que nunca es necesario rep.
Angs
Puede usar x<-0:10/10;f<-function(x){-2};10^-2*sapply(x,f)si f (x) no está garantizado para ser vectorizado o simplemente 10^-2*f(x)si festá vectorizado ( laplace(Vectorize(function(x)2),0,0,10)
Angs
No use eval, dar fcomo una función adecuada.
Angs
2

Haskell, 195 168 bytes

import Numeric.LinearAlgebra
p f s e n|m<-[0..]!!n=((n><n)(([1,0]:([3..n]>>[[-1,2,-1]])++[[0,1]])>>=(++(0<$[3..n]))))<\>(col$s:map((/(m-1)^2).f.(/(m-1)))[1..m-2]++[e])

La legibilidad recibió un gran golpe. Sin golf:

laplace f start end _N = linearSolveLS _M y
  where
  n = fromIntegral _N
  _M = (_N><_N) --construct n×n matrix from list
        ( ( [1,0]           --make a list of [1,0]
          : ([3.._N]>>[[-1,2,-1]]) --         (n-2)×[-1,2,-1]
          ++ [[0,1]])       --               [0,1]
        >>= (++(0<$[3.._N])) --append (n-2) zeroes and concat
        )                   --(m><n) discards the extra zeroes at the end
  h  = 1/(n-1) :: Double
  y  = asColumn . fromList $ start : map ((*h^2).f.(*h)) [1..n-2] ++ [end]

TODO: Imprimiendo en 83 71 bytes.

Déjame ver:

import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo

D'oh!

Angs
fuente
No sé mucho sobre Haskell, pero quizás la solución sin álgebra matricial podría ser más corta, agregué una segunda implementación de muestra.
Karl Napf
@KarlNapf no se acerca mucho Esto es solo semi-golf pero tiene que usar muchas funciones detalladas integradas. Con álgebra matricial, la mayor parte del código está construyendo la matriz (64 bytes) y la importación (29 bytes). El residual y el jacobi ocupan mucho espacio.
Angs
Bueno, muy mal, pero valió la pena intentarlo.
Karl Napf
1

Axioma, 579 460 bytes

l(w,y)==(r:=0;for i in 1..y|index?(i,w)repeat r:=i;r)
g(z:EQ EXPR INT,y:BasicOperator,a0:Float,a1:Float,a2:Float):Float==(r:=digits();digits(r+30);q:=seriesSolve(z,y,x=0,[a,b])::UTS(EXPR INT,x,0);w:=eval(q,0);s:=l(w,r+30);o:=solve([w.s=a0,eval(q,1).s=a1]::List(EQ POLY Float),[a,b]);v:=eval(eval(eval(q,a2).s,o.1.1),o.1.2);digits(r);v)
m(z:EXPR INT,a0:Float,a1:Float,n:INT):List Float==(n:=n-1;y:=operator 'y;r:=[g(D(y x,x,2)=-z,y,a0,a1,i/n)for i in 0..n];r)

deshazte de él y prueba

Len(w,y)==(r:=0;for i in 1..y|index?(i,w)repeat r:=i;r)

-- g(z,a0,a1,a2)
-- Numeric solve z as f(y''(x),y'(x),y(x))=g(x) with ini conditions y(0)=a0   y(1)=a1 in x=a2
NSolve2order(z:EQ EXPR INT,y:BasicOperator,a0:Float,a1:Float,a2:Float):Float==
      r:=digits();digits(r+30)
      q:=seriesSolve(z,y,x=0,[a,b])::UTS(EXPR INT,x,0)
      w:=eval(q,0);s:=Len(w,r+30)
      o:=solve([w.s=a0,eval(q,1).s=a1]::List(EQ POLY Float),[a,b])
      v:=eval(eval(eval(q,a2).s,o.1.1),o.1.2);digits(r)
      v

InNpoints(z:EXPR INT,a0:Float,a1:Float,n:INT):List Float==(n:=n-1;y:=operator 'y;r:=[NSolve2order(D(y x,x,2)=-z,y,a0,a1,i/n)for i in 0..n];r)

la función para la pregunta es m (,,,) el código anterior se coloca en el archivo "file.input" y se carga en Axiom. El resultado depende de la función digits ().

si alguien piensa que no es golf => él o ella puede mostrar cómo hacerlo ... gracias

PD

Parecen los 6 números después del. para e ^ (x ^ 2) no están bien aquí o en los ejemplos, pero aquí aumento los dígitos pero los números no cambian ... para mí significa que los números en el ejemplo están mal. ¿Por qué todos los demás no han mostrado sus números?

para sin (2 *% pi * x) también hay problemas

"Aquí la solución exacta es u = (sin (2 * π * x)) / (4 * π ^ 2) +1" copié la solución exacta para x = 1/19:

              (sin(2*π/19))/(4*π^2)+1

en WolframAlpha https://www.wolframalpha.com/input/?i=(sin(2% CF% 80% 2F19))% 2F (4 % CF% 80% 5E2)% 2B1 resulta

1.008224733636964333380661957738992274267070440829381577926...
1.0083001
  1234
1.00822473

1.0083001 propuesto como resultado es diferente en el cuarto dígito del resultado real 1.00822473 ... (y no sexto)

-- in interactive mode
(?) -> )read  file
(10) -> digits(9)
   (10)  10
                                                        Type: PositiveInteger
(11) -> m(-2,0,0,10)
   (11)
   [0.0, - 0.0987654321, - 0.172839506, - 0.222222222, - 0.24691358,
    - 0.24691358, - 0.222222222, - 0.172839506, - 0.098765432, 0.0]
                                                             Type: List Float
(12) -> m(10*x,0,1,15)
   (12)
   [0.0, 0.189868805, 0.376093294, 0.555029155, 0.72303207, 0.876457726,
    1.01166181, 1.125, 1.21282799, 1.27150146, 1.29737609, 1.28680758,
    1.2361516, 1.14176385, 1.0]
                                                             Type: List Float
(13) -> m(sin(2*%pi*x),1,1,20)
   (13)
   [1.0, 1.00822473, 1.01555819, 1.02120567, 1.0245552, 1.02524378, 1.02319681,
    1.0186361, 1.01205589, 1.00416923, 0.99583077, 0.987944112, 0.981363896,
    0.976803191, 0.97475622, 0.975444804, 0.978794326, 0.98444181, 0.991775266,
    1.0]
                                                         Type: List Float
(14) -> m(exp(x^2),0,0,30)
   (14)
   [0.0, 0.0202160702, 0.0392414284, 0.0570718181, 0.0737001105, 0.0891162547,
    0.103307204, 0.116256821, 0.127945761, 0.138351328, 0.147447305,
    0.155203757, 0.161586801, 0.166558343, 0.170075777, 0.172091643,
    0.172553238, 0.171402177, 0.168573899, 0.163997099, 0.157593103,
    0.149275146, 0.13894757, 0.126504908, 0.111830857, 0.0947971117,
    0.0752620441, 0.0530692118, 0.0280456602, - 0.293873588 E -38]
                                                             Type: List Float
RosLuP
fuente
La solución numérica difiere de la solución exacta porque el FDM aquí es solo de segundo orden, lo que significa que solo los polinomios hasta el orden 2 pueden representarse exactamente. Entonces, solo el f=-2ejemplo tiene una solución analítica y numérica coincidente.
Karl Napf
aquí la solución numérica parece correcta si cambio los dígitos a 80 o 70 -> g (sin (2 *% pi * x), 1,1,1 / 19) 1.0082247336 3696433338 0661957738 9922742670 7044082938 1577926908 950765832 el otro número 1.0082247336 3696433338 0667426773899 7044082938 1577926 ...
RosLuP