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!

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]

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

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:

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)

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!

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