Rozwiąż równanie Laplace'a

13

Wprowadzenie do matematyki numerycznej

To jest „Witaj, świecie!” PDE (równania różniczkowe cząstkowe). Równanie Laplace'a lub dyfuzji pojawia się często w fizyce, na przykład równanie ciepła, deformacja, dynamika płynów itp. Tak jak w prawdziwym życiu jest 3D, ale chcemy powiedzieć „Cześć, świecie!” i nie zaśpiewajcie „99 butelek piwa ...” to zadanie podano w 1D. Możesz to zinterpretować jako gumową szatę przywiązaną do ściany na obu końcach z użyciem pewnej siły.

W [0,1]domenie znajdź funkcję udla danej funkcji źródłowej fi wartości brzegowych u_Li u_Rtaką, aby:

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

u'' oznacza drugą pochodną u

Można to rozwiązać czysto teoretycznie, ale Twoim zadaniem jest rozwiązanie go liczbowo na dyskretnej domenie x dla Npunktów:

  • x = {i/(N-1) | i=0..N-1}lub w oparciu o 1:{(i-1)/(N-1) | i=1..N}
  • h = 1/(N-1) to odstępy

Wejście

  • f jako funkcja lub wyrażenie lub ciąg
  • u_L, u_Rjako wartości zmiennoprzecinkowe
  • N jako liczba całkowita> = 2

Wynik

  • Tablica, lista, jakiś oddzielny ciąg utakich znakówu_i == u(x_i)

Przykłady

Przykład 1

Wejście: f = -2, u_L = u_R = 0, N = 10(nie ma f=-2źle, to nie jest wartością stałą, ale funkcja, która powraca -2do wszystkich xTo jak stałą grawitacji życie naszej liny.).

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

Istnieje łatwe dokładne rozwiązanie: u = -x*(1-x)

Przykład 2

Wejście: f = 10*x, u_L = 0 u_R = 1, N = 15(Tu jest dużo wiatr z prawej strony)

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

Dokładne rozwiązanie tego stanu: u = 1/3*(8*x-5*x^3)

Przykład 3

Wejście: f = sin(2*pi*x), u_L = u_R = 1, N = 20(Ktoś złamał grawitacji lub nie jest rodzajem górę iz wiatrem)

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

Oto dokładne rozwiązanie u = (sin(2*π*x))/(4*π^2)+1

Przykład 4

Wejście: f = exp(x^2), u_L = u_R = 0,N=30

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

Zwróć uwagę na niewielką asymetrię

FDM

Jedną z możliwych metod rozwiązania tego problemu jest metoda różnic skończonych :

  • przepisz -u_i'' = f_ijako
  • (-u_{i-1} + 2u_i - u{i+1})/h² = f_i co jest równe
  • -u_{i-1} + 2u_i - u{i+1} = h²f_i
  • Skonfiguruj równania:

  • Które są równe równaniu macierz-wektor:

  • Rozwiąż to równanie i wyślij u_i

Jedna implementacja tego do demonstracji w Pythonie:

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)

Alternatywne wdrożenie bez Matrycy Algebry (przy użyciu metody 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

Możesz jednak użyć dowolnej innej metody, aby rozwiązać równanie Laplace'a. Jeśli używasz metody iteracyjnej, powinieneś iterować do reszty |b-Au|<1e-6, bbędąc wektorem po prawej stronieu_L,f_1h²,f_2h²,...

Notatki

W zależności od metody rozwiązania przykłady mogą nie zostać rozwiązane dokładnie dla podanych rozwiązań. Przynajmniej N->infinitybłąd powinien zbliżyć się do zera.

Standardowe luki są niedozwolone , dozwolone są wbudowane PDE.

Premia

Premia w wysokości -30% za wyświetlenie rozwiązania graficznego lub ASCII-art.

Zwycięski

To jest codegolf, więc wygrywa najkrótszy kod w bajtach!

Karl Napf
źródło
Polecam dodanie przykładu, który nie jest analitycznie możliwy do rozwiązania, np f(x) = exp(x^2). Z.
flawr
@flawr Jasne, ma rozwiązanie, jednak występuje funkcja błędu.
Karl Napf,
1
Przepraszam, to może być niewłaściwe wyrażenie, czy „niepierwiastkowa funkcja pierwotna” może być bardziej odpowiednia? Mam na myśli funkcje takie jak log(log(x))lub, sqrt(1-x^4)które mają całkę, która jednak nie jest wyrażalna w funkcjach elementarnych.
flawr
@ flawr Nie jest w porządku, funkcja błędu nie jest elementarna, chciałem tylko powiedzieć, że istnieje sposób analitycznego wyrażenia rozwiązania, ale u(x) = 1/2 (-sqrt(π) x erfi(x)+sqrt(π) erfi(1) x+e^(x^2)-e x+x-1)nie jest dokładnie obliczalny.
Karl Napf,
Dlaczego iterować do 1e-6, a nie iterować do 1e-30?
RosLuP,

Odpowiedzi:

4

Mathematica, 52,5 bajtów (= 75 * (1-30%))

+0,7 bajtów na komentarz @flawr.

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

To drukuje dane wyjściowe.

na przykład

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

wprowadź opis zdjęcia tutaj

Wyjaśnienie

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

Rozwiąż dla funkcji u.

Subdivide@#4

Subdivide przedział [0,1] na części N (4. wejście).

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

Odwzoruj una wynik Subdivide.

ListPlot[ ... ]

Wykreśl wynik końcowy.

Rozwiązanie inne niż graficzne: 58 bajtów

u/@Subdivide@#4/.NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]&
JungHwan Min
źródło
To nie działa w przypadkuf(x) = exp(x^2)
flawr
Być może możesz chcieć użyć NDSolvedo ogólnego przypadku rozwiązań nie elementarnych.
flawr
6

Matlab, 84, 81,2 79,1 bajtów = 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)

Zauważ, że w tym przykładzie używam wektorów wierszy, co oznacza, że ​​macierz Ajest transponowana. fjest traktowany jako uchwyt funkcyjny, a,bsą to ograniczenia Dirichleta po lewej / prawej stronie.

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

Na przykład f = 10*x, u_L = 0 u_R = 1, N = 15powoduje to:

wada
źródło
3

R, 123,2 102,9 98,7 bajtów (141-30%)

Edycja: Zaoszczędź garść bajtów dzięki @Angs!

Jeśli ktoś chce edytować zdjęcia, możesz to zrobić. Jest to w zasadzie adaptacja R zarówno opublikowanych wersji Matlaba, jak i Pythona.

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))}

Niegolfowane i wyjaśnione:

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 
}

Przykłady i przypadki testowe:

Nazwaną i nie golfową funkcję można wywołać za pomocą:

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)

Zauważ, że fargument jest funkcją R.

Aby uruchomić wersję golfową, wystarczy użyć (function(...){...})(args)

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

Billywob
źródło
Myślę, że możesz zrezygnować z is.numeric(f)czeku, jeśli zadeklarujesz fjako funkcję, nie jest wymagane przekazywanie go bezpośrednio w wywołaniu funkcji do solvera.
Karl Napf,
Ach, rozumiem, nie wiedziałem, że jest różnica między tymi dwoma. Cóż, jeśli jest krótszy, możesz zmodyfikować swój solver tak, aby akceptował fjako funkcję, więc nie musisz sprawdzać, czy jest to stała (funkcja).
Karl Napf,
1
@Billywob, nigdy nie trzeba fbyć cyfrowym . f = (function(x)-2)działa dla pierwszego przykładu, więc nigdy nie ma takiej potrzeby rep.
Angs
Można użyć x<-0:10/10;f<-function(x){-2};10^-2*sapply(x,f), jeśli f (x) nie jest quaranteed być wektorowy lub po prostu 10^-2*f(x)jeśli fjest wektorowy ( laplace(Vectorize(function(x)2),0,0,10)
Angs
Nie używaj eval, daj fjako właściwą funkcję.
Angs,
2

Haskell, 195 168 bajtów

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

Czytelność spotkała się z dużym hitem. Nie golfowany:

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]

DO ZROBIENIA: Drukowanie w 83 71 bajtach.

Daj mi zobaczyć:

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

Nie!

Angs
źródło
Nie wiem dużo o Haskell, ale być może rozwiązanie bez algebry macierzy może być krótsze, dodałem drugą przykładową implementację.
Karl Napf,
@KarlNapf nie podchodzi bardzo blisko Jest to tylko gra w golfa częściowo, ale musi ona używać wielu pełnych wbudowanych funkcji. W przypadku algebry macierzowej większość kodu buduje macierz (64 bajty) i import (29 bajtów). Resztki i Jacobi zajmują sporo miejsca.
Angs,
Cóż, szkoda, ale warto było spróbować.
Karl Napf,
1

Axiom, 579 460 bajtów

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)

rozkop go i przetestuj

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)

funkcją pytania jest m (,,,) powyższy kod jest umieszczony w pliku „file.input” i załadowany do Axiom. Wynik zależy od funkcji digits ().

jeśli ktoś myśli, że to nie gra w golfa => on lub ona może pokazać, jak to zrobić ... dzięki

PS

wydaje się, że 6 liczb po. dla e ^ (x ^ 2) nie są tutaj ani w przykładach, ale tutaj zwiększam cyfry, ale liczby się nie zmieniają ... dla mnie oznacza to, że liczby w przykładzie są nieprawidłowe. Dlaczego wszyscy inni nie pokazali swoich liczb?

w przypadku sin (2 *% pi * x) występują również problemy

„Tutaj dokładne rozwiązanie to u = (sin (2 * π * x)) / (4 * π ^ 2) +1” skopiowałem dokładne rozwiązanie dla x = 1/19:

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

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

1.008224733636964333380661957738992274267070440829381577926...
1.0083001
  1234
1.00822473

1.0083001 zaproponowany jako wynik różni się czwartą cyfrą od rzeczywistego wyniku 1.00822473 ... (a nie 6.)

-- 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
źródło
Rozwiązanie numeryczne różni się od rozwiązania dokładnego, ponieważ FDM jest tutaj tylko drugim rzędem, co oznacza, że ​​tylko wielomiany do rzędu 2 mogą być dokładnie reprezentowane. Tak więc tylko f=-2przykład ma pasujące rozwiązanie analityczne i numeryczne.
Karl Napf,
tutaj rozwiązanie numeryczne wydaje się być prawidłowe, jeśli zmienię cyfry na 80 lub 70 -> g (sin (2 *% pi * x), 1,1,1 / 19) 1.0082247336 3696433338 0661957738 9922742670 7044082938 1577926908 950765832 inny numer 1.0082247336 36964333333 99 0 7044082938 1577926 ...
RosLuP