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!

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]

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]&

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:

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)

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!

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