Resolva a equação de Laplace

13

Introdução à Matemática Numérica

Este é o "Olá, mundo!" de PDEs (Equações Diferenciais Parciais). A equação de Laplace ou difusão aparece frequentemente em física, por exemplo, equação de calor, deformação, dinâmica de fluidos, etc ... Como a vida real é 3D, mas queremos dizer "Olá, mundo!" e não cantar "99 garrafas de cerveja, ...", esta tarefa é dada em 1D. Você pode interpretar isso como uma túnica de borracha amarrada a uma parede nas duas extremidades com alguma força aplicada a ela.

Em um [0,1]domínio, encontre uma função upara determinada função de origem fe valores de limite u_Le de u_Rforma que:

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

u'' denota a segunda derivada de u

Isso pode ser resolvido puramente teórico, mas sua tarefa é resolvê-lo numericamente em um domínio discretizado x para obter Npontos:

  • x = {i/(N-1) | i=0..N-1}ou com base em 1:{(i-1)/(N-1) | i=1..N}
  • h = 1/(N-1) é o espaçamento

Entrada

  • f como função ou expressão ou string
  • u_L, u_Rcomo valores de ponto flutuante
  • N como inteiro> = 2

Resultado

  • Array, List, algum tipo de string separada de utal modo queu_i == u(x_i)

Exemplos

Exemplo 1

Entrada: f = -2, u_L = u_R = 0, N = 10(não pegue f=-2errado, não é um valor, mas uma função constante que os retornos -2para todo xÉ como uma força de gravidade constante em nosso corda..)

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

Existe uma solução exata e fácil: u = -x*(1-x)

Exemplo 2

Entrada: f = 10*x, u_L = 0 u_R = 1, N = 15(Aqui há um monte de contra o vento no lado direito)

Resultado: [ 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.]

A solução exata para esses estados: u = 1/3*(8*x-5*x^3)

Exemplo 3

Entrada: f = sin(2*pi*x), u_L = u_R = 1, N = 20(Alguém quebrou gravidade ou há uma espécie de cima e na direção do vento)

Resultado: [ 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.]

Aqui a solução exata é u = (sin(2*π*x))/(4*π^2)+1

Exemplo 4

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

Resultado: [ 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. ]

Observe a leve simetria

FDM

Um método possível para resolver isso é o método das diferenças finitas :

  • reescrever -u_i'' = f_icomo
  • (-u_{i-1} + 2u_i - u{i+1})/h² = f_i que é igual
  • -u_{i-1} + 2u_i - u{i+1} = h²f_i
  • Configure as equações:

  • Que são iguais a uma equação de vetor de matriz:

  • Resolva esta equação e produza o u_i

Uma implementação disso para demonstração em 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)

Implementação alternativa sem Álgebra Matricial (usando o 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

No entanto, você pode usar qualquer outro método para resolver a equação de Laplace. Se você usar um método iterativo, iterará até o residual |b-Au|<1e-6, bsendo o vetor do lado direitou_L,f_1h²,f_2h²,...

Notas

Dependendo do seu método de solução, você pode não resolver os exemplos exatamente para as soluções fornecidas. Pelo menos para N->infinityo erro deve se aproximar de zero.

As brechas padrão são proibidas , os internos para PDEs são permitidos.

Bônus

Um bônus de -30% para exibir a solução, gráfica ou ASCII-art.

Ganhando

Isso é codegolf, então o código mais curto em bytes vence!

code-golf math number Karl Napf
fonte
Eu recomendo adicionar um exemplo que não seja analiticamente solucionável, por exemplo, com f(x) = exp(x^2).
flawr
@ flawr Claro, ele tem uma solução, no entanto, a função de erro está envolvida.
194 Karl Napf
1
Desculpe, talvez essa tenha sido a expressão errada. O "antiderivativo não elementar" pode ser mais adequado? Quero dizer funções como log(log(x))ou sqrt(1-x^4)que possuem uma integral, que, no entanto, não é expressável em funções elementares.
flawr
@ flawr Não, tudo bem, a função de erro não é elementar, eu só queria dizer que há uma maneira de expressar a solução analiticamente, mas u(x) = 1/2 (-sqrt(π) x erfi(x)+sqrt(π) erfi(1) x+e^(x^2)-e x+x-1)não é exatamente calculável.
194 Karl Napf
Por que iterar até 1e-6 e não iterar até 1e-30?
RosLuP

Respostas:

4

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

+0,7 bytes por comentário de @flawr.

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

Isso plota a saída.

por exemplo

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

insira a descrição da imagem aqui

Explicação

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

Resolva para a função u.

Subdivide@#4

Subdivide o intervalo [0,1] em N (quarta entrada) partes.

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

Mapeie upara a saída de Subdivide.

ListPlot[ ... ]

Traçar o resultado final.

Solução não gráfica: 58 bytes

u/@Subdivide@#4/.NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]&
JungHwan Min
fonte
Isso não funciona paraf(x) = exp(x^2)
flawr 7/16
Talvez você queira usar NDSolveno caso geral de soluções não elementares.
flawr
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)

Observe que, neste exemplo, eu uso vetores de linha, isso significa que a matriz Aé transposta. fé tomado como identificador de função, a,bsão as restrições de Dirichlet do lado esquerdo / direito.

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 o exemplo, f = 10*x, u_L = 0 u_R = 1, N = 15isso resulta em:

flawr
fonte
3

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

Editar: salvou um punhado de bytes graças a @Angs!

Se alguém quiser editar as imagens, sinta-se à vontade para fazê-lo. Esta é basicamente uma adaptação em R das versões matlab e python postadas.

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 & explicou:

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 
}

Exemplos e casos de teste:

A função nomeada e ungolfed pode ser chamada 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)

Observe que o fargumento é uma função R.

Para executar a versão golfed simplesmente use (function(...){...})(args)

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

Billywob
fonte
Eu acho que você pode cancelar a is.numeric(f)verificação se você declarar fcomo função, não é necessário passá-lo diretamente na chamada de função para o solucionador.
194 Karl Napf
Ah entendo, eu não sabia que havia uma diferença entre os dois. Bem, se ele for mais curto, você poderá modificar seu solucionador para aceitar fcomo uma função, para não precisar verificar se o caso é uma constante (função).
194 Karl Napf
1
@ Billywob, não há necessidade de fser numérico. f = (function(x)-2)funciona para o primeiro exemplo, por isso nunca é necessário rep.
Angs
Você pode usar x<-0:10/10;f<-function(x){-2};10^-2*sapply(x,f), se f (x) não é quaranteed a ser vetorizado ou apenas 10^-2*f(x)se festá vetorizado ( laplace(Vectorize(function(x)2),0,0,10)
Angs
Não use eval, dê fcomo uma função adequada.
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])

A legibilidade levou bastante sucesso. Ungolfed:

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: Imprimindo em 83 71 bytes.

Deixa eu ver:

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

D'oh!

Angs
fonte
Não sei muito sobre Haskell, mas talvez a solução sem álgebra matricial possa ser mais curta, adicionei uma segunda implementação de amostra.
9136 Karl Napf
@KarlNapf não chega muito perto Isso é apenas semi-golfe, mas tem que usar muitas funções internas detalhadas. Com a álgebra matricial, a maior parte do código está construindo a matriz (64 bytes) e a importação (29 bytes). O residual e o jacobi ocupam bastante espaço.
Angs
Bem, muito ruim, mas valeu a pena tentar.
Karl Napf 7/11
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)

ungolf-lo e testar

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)

a função da pergunta é m (,,,) o código acima é colocado no arquivo "file.input" e carregado no Axiom. O resultado depende da função digits ().

se alguém achar que não é jogador de golfe => ele ou ela pode mostrar como fazer isso ... obrigado

PS

parece que os 6 números depois da. para e ^ (x ^ 2) não está ok aqui ou nos exemplos, mas aqui eu aumento os dígitos, mas os números não mudam ... para mim, isso significa que os números no exemplo estão errados. Por que todos os outros não mostraram seus números?

para o pecado (2 *% pi * x) também existem problemas

"Aqui a solução exata é u = (sin (2 * π * x)) / (4 * π ^ 2) +1" copiei a solução exata para x = 1/19: 

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

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

1.008224733636964333380661957738992274267070440829381577926...
1.0083001
  1234
1.00822473

1.0083001 proposto como resultado é diferente no quarto dígito do resultado real 1.00822473 ... (e não 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
fonte
A solução numérica difere da solução exata porque o FDM aqui é apenas de segunda ordem, o que significa que apenas polinômios até a ordem 2 podem ser representados exatamente. Portanto, apenas o f=-2exemplo tem uma solução analítica e numérica correspondente.
11136 Karl Napf
aqui a solução numérica parece ok se eu alterar os dígitos para 80 ou 70 -> g (sin (2 *% pi * x), 1,1,1 / 19) 7044082938 1577926 ...
RosLuP 13/11