Lösen Sie die Laplace-Gleichung

13

Einführung in die Numerische Mathematik

Dies ist das "Hallo, Welt!" von PDEs (Partielle Differentialgleichungen). Die Laplace- oder Diffusionsgleichung wird in der Physik häufig verwendet, z. B. Wärmegleichung, Verformung, Fluiddynamik usw. Da das reale Leben 3D ist, möchten wir jedoch "Hallo, Welt!" und nicht singen "99 flaschen bier, ..." diese aufgabe ist in 1D gegeben. Sie können dies als eine Gummimantel interpretieren, die an beiden Enden mit etwas Kraft an einer Wand befestigt ist.

[0,1]Suchen Sie in einer Domain eine Funktion ufür bestimmte Quellfunktionen fund Grenzwerte, u_Lund zwar u_Rso, dass:

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

u'' bezeichnet die zweite Ableitung von u

Dies kann rein theoretisch gelöst werden, aber Ihre Aufgabe ist es, es numerisch in einer diskretisierten Domäne x nach NPunkten zu lösen :

  • x = {i/(N-1) | i=0..N-1}oder 1-basiert:{(i-1)/(N-1) | i=1..N}
  • h = 1/(N-1) ist der Abstand

Eingang

  • f als Funktion oder Ausdruck oder Zeichenfolge
  • u_L, u_RAls Gleitkommawerte
  • N als ganze Zahl> = 2

Ausgabe

  • Array, List, eine Art separater String udavonu_i == u(x_i)

Beispiele

Beispiel 1

Input: f = -2, u_L = u_R = 0, N = 10(nehmen Sie nicht f=-2falsch, es ist kein Wert , sondern eine konstante Funktion , dass die Renditen -2für alle xEs ist wie eine konstante Schwerkraft auf unserem Seil..)

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

Es gibt eine einfache exakte Lösung: u = -x*(1-x)

Beispiel 2

Input: f = 10*x, u_L = 0 u_R = 1, N = 15(Hier gibt es eine Menge gegen den Wind auf der rechten Seite)

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

Die genaue Lösung hierfür lautet: u = 1/3*(8*x-5*x^3)

Beispiel 3

Input: f = sin(2*pi*x), u_L = u_R = 1, N = 20(jemand brach die Schwerkraft , oder es ist eine Art von Up- und Abwind)

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

Hier ist die genaue Lösung u = (sin(2*π*x))/(4*π^2)+1

Beispiel 4

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

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

Beachten Sie die leichte Unsymmetrie

FDM

Eine mögliche Methode, um dies zu lösen, ist die Finite-Differenz-Methode :

  • umschreiben -u_i'' = f_ials
  • (-u_{i-1} + 2u_i - u{i+1})/h² = f_i was gleich ist
  • -u_{i-1} + 2u_i - u{i+1} = h²f_i
  • Stellen Sie die Gleichungen auf:

  • Welche sind gleich einer Matrix-Vektor-Gleichung:

  • Lösen Sie diese Gleichung und geben Sie die aus u_i

Eine Implementierung davon zur Demonstration in 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)

Alternative Implementierung ohne Matrixalgebra (unter Verwendung der Jacobi-Methode )

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

Sie können jedoch eine andere Methode verwenden, um die Laplace-Gleichung zu lösen. Wenn Sie eine iterative Methode verwenden, sollten Sie bis zum Residuum iterieren |b-Au|<1e-6, wobei bdies der Vektor auf der rechten Seite istu_L,f_1h²,f_2h²,...

Anmerkungen

Abhängig von Ihrer Lösungsmethode lösen Sie die Beispiele möglicherweise nicht genau mit den angegebenen Lösungen. Zumindest für N->infinityden Fehler sollte sich Null nähern.

Standard-Schlupflöcher sind nicht zulässig, eingebaute PDEs sind zulässig.

Bonus

Ein Bonus von -30% für die Anzeige der Lösung, entweder grafisch oder ASCII-artig.

Gewinnen

Das ist Codegolf, also gewinnt der kürzeste Code in Bytes!

code-golf math number Karl Napf
quelle
Ich empfehle, ein Beispiel hinzuzufügen, das analytisch nicht lösbar ist, zB mit f(x) = exp(x^2).
Fehler
@flawr Klar, es gibt eine Lösung, jedoch ist die Fehlerfunktion involviert.
Karl Napf
1
Sorry, das war vielleicht der falsche Ausdruck, könnte "nicht-elementares Antiderivativ" besser geeignet sein? Ich meine Funktionen wie log(log(x))odersqrt(1-x^4) die ein Integral haben, das sich aber in Elementarfunktionen nicht ausdrücken lässt.
Fehler
@flawr Nein, es ist in Ordnung, die Fehlerfunktion ist nicht elementar, ich wollte nur sagen, es gibt eine Möglichkeit, die Lösung aber analytisch auszudrücken u(x) = 1/2 (-sqrt(π) x erfi(x)+sqrt(π) erfi(1) x+e^(x^2)-e x+x-1) nicht genau berechenbar ist.
Karl Napf
Warum bis 1e-6 iterieren und nicht bis 1e-30 iterieren?
RosLuP

Antworten:

4

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

+0,7 Bytes pro Kommentar von @flawr.

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

Dies zeichnet die Ausgabe.

z.B

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

Bildbeschreibung hier eingeben

Erläuterung

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

Löse für die Funktion u.

Subdivide@#4

Subdivide das Intervall [0,1] in N (4. Eingabe) Teile.

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

Map uauf die Ausgabe von Subdivide.

ListPlot[ ... ]

Zeichnen Sie das Endergebnis.

Lösung ohne Grafik: 58 Byte

u/@Subdivide@#4/.NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]&
JungHwan min
quelle
Dies funktioniert nicht fürf(x) = exp(x^2)
Fehler
Vielleicht möchten Sie NDSolvefür den allgemeinen Fall nicht-elementarer Lösungen verwenden.
Fehler
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)

Beachten Sie, dass in diesem Beispiel ich Zeilenvektoren verwende, das bedeutet, dass die Matrix Atransponiert ist. fwird als Funktionshandle genommen,a,b sind die Dirichlet-Beschränkungen links / rechts.

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

Für das Beispiel f = 10*x, u_L = 0 u_R = 1, N = 15ergibt dies:

Fehler
quelle
3

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

Bearbeiten: Eine Handvoll Bytes dank @Angs gespeichert!

Wenn jemand die Bilder bearbeiten möchte, kann er dies gerne tun. Dies ist im Grunde eine R-Anpassung sowohl der Matlab- als auch der Python-Version.

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 & erklärt:

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 
}

Beispiel & Testfälle:

Die Funktion named und ungolfed kann folgendermaßen aufgerufen werden:

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)

Notiere dass der f Argument eine R-Funktion ist.

Um die Golfversion zu starten, verwenden Sie einfach (function(...){...})(args)

Bildbeschreibung hier eingeben Bildbeschreibung hier eingeben Bildbeschreibung hier eingeben Bildbeschreibung hier eingeben

Billywob
quelle
Ich denke, Sie können die is.numeric(f)Prüfung verwerfen, wenn Sie fals Funktion deklarieren . Es ist nicht erforderlich, sie direkt im Funktionsaufruf an den Solver zu übergeben.
Karl Napf
Ah ich verstehe, ich wusste nicht, dass es einen Unterschied zwischen diesen beiden gibt. Wenn es kürzer ist, können Sie Ihren Solver fso ändern, dass er als Funktion akzeptiert wird, sodass Sie nicht prüfen müssen, ob es sich um eine Konstante (Funktion) handelt.
Karl Napf
1
@ Billywob Es ist nicht nötig f, jemals numerisch zu sein. f = (function(x)-2)funktioniert für das erste Beispiel, es ist also nie nötig rep.
Angs
Sie können verwenden, x<-0:10/10;f<-function(x){-2};10^-2*sapply(x,f)wenn f (x) nicht als vektorisiert unter Quarantäne gestellt wird, oder nur, 10^-2*f(x)wenn fvektorisiert ist ( laplace(Vectorize(function(x)2),0,0,10)
Angs
Verwenden Sie eval nicht, sondern geben Sie es fals richtige Funktion an.
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])

Die Lesbarkeit hat einiges gekostet. 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: Drucken in 83 71 Bytes.

Lass mich sehen:

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

D'oh!

Angs
quelle
Ich weiß nicht viel über Haskell, aber vielleicht ist die Lösung ohne Matrixalgebra kürzer. Ich habe eine zweite Beispielimplementierung hinzugefügt.
Karl Napf
@KarlNapf kommt nicht sehr nahe Dies ist nur halbgolfen, aber es muss viele ausführliche integrierte Funktionen verwenden. Bei der Matrixalgebra besteht der größte Teil des Codes aus der Matrix (64 Byte) und dem Import (29 Byte). Das Residuum und das Jacobi nehmen ziemlich viel Platz ein.
Angs
Schade, aber es war einen Versuch wert.
Karl Napf
1

Axiom, 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)

entgolfen und testen

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)

die funktion für die frage ist m (,,,) der obige code wird in die datei "file.input" gestellt und in axiom geladen. Das Ergebnis hängt von der Funktion digits () ab.

Wenn jemand denkt, dass es kein Golfspiel ist => kann er oder sie zeigen, wie es geht ... danke

PS

es scheint die 6 Zahlen nach dem. für e ^ (x ^ 2) sind hier oder in den Beispielen nicht in Ordnung, aber hier erhöhe ich die Ziffern, aber die Zahlen ändern sich nicht ... für mich bedeutet das, dass die Zahlen im Beispiel falsch sind. Warum haben alle anderen ihre Nummern nicht gezeigt?

für sin (2 *% pi * x) gibt es auch probleme

"Hier ist die exakte Lösung u = (sin (2 * π * x)) / (4 * π ^ 2) +1" Ich habe die exakte Lösung für x = 1/19 kopiert: 

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

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

1.008224733636964333380661957738992274267070440829381577926...
1.0083001
  1234
1.00822473

1.0083001 als Ergebnis vorgeschlagen unterscheidet sich in der 4. Ziffer vom tatsächlichen Ergebnis 1.00822473 ... (und nicht 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
quelle
Die numerische Lösung unterscheidet sich von der exakten Lösung, da die FDM hier nur zweiter Ordnung ist, dh nur Polynome bis zur Ordnung 2 exakt dargestellt werden können. So hat nur das f=-2Beispiel eine passende analytische und numerische Lösung.
Karl Napf
Hier scheint die numerische Lösung in Ordnung zu sein, wenn ich die Ziffern auf 80 oder 70 ändere -> g (sin (2 *% pi * x), 1,1,1 / 19) 7044082938 1577926 ...
RosLuP