Pengantar Matematika Numerik

Ini adalah "Halo, Dunia!" PDE (Persamaan Diferensial Parsial). Persamaan Laplace atau Difusi sering muncul dalam Fisika, misalnya Persamaan Panas, Deforming, Fluid Dynamics, dll ... Karena kehidupan nyata adalah 3D tetapi kami ingin mengatakan "Halo, Dunia!" dan tidak menyanyikan "99 botol bir, ..." tugas ini diberikan dalam 1D. Anda dapat menafsirkan ini sebagai jubah karet yang diikat ke dinding di kedua ujungnya dengan beberapa kekuatan yang diterapkan padanya.

Pada [0,1]domain, temukan fungsi uuntuk fungsi sumber fdan nilai batas yang diberikan u_Ldan u_Rsedemikian rupa sehingga:

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

u'' menunjukkan turunan kedua dari u

Ini dapat dipecahkan sepenuhnya teoretis tetapi tugas Anda adalah menyelesaikannya secara numerik pada domain yang didiskritisasi x untuk Npoin:

• x = {i/(N-1) | i=0..N-1}atau berbasis 1:{(i-1)/(N-1) | i=1..N}
• h = 1/(N-1) adalah spasi

Memasukkan

• f sebagai fungsi atau ekspresi atau string
• u_L, u_Rsebagai nilai floating point
• N sebagai integer> = 2

Keluaran

• Array, List, semacam string terpisah usehinggau_i == u(x_i)

Contohnya

Contoh 1

Input: f = -2, u_L = u_R = 0, N = 10(Do tidak mengambil f=-2salah, itu bukan nilai tapi fungsi konstan yang kembali -2untuk semua xHal ini seperti gaya gravitasi konstan pada tali kami..)

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

Ada solusi tepat yang mudah: u = -x*(1-x)

Contoh 2

Input: f = 10*x, u_L = 0 u_R = 1, N = 15(Di sini ada banyak melawan angin di sisi kanan)

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

Solusi tepat untuk kondisi ini: u = 1/3*(8*x-5*x^3)

Contoh 3

Input: f = sin(2*pi*x), u_L = u_R = 1, N = 20(Seseorang masuk gravitasi atau ada semacam atas dan melawan arah angin)

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

Di sini solusinya adalah u = (sin(2*π*x))/(4*π^2)+1

Contoh 4

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

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

Perhatikan sedikit ketidaksimetrisan

FDM

Salah satu metode yang mungkin untuk menyelesaikan ini adalah Metode Perbedaan Hingga :

• tulis ulang -u_i'' = f_isebagai
• (-u_{i-1} + 2u_i - u{i+1})/h² = f_i yang sama dengan
• -u_{i-1} + 2u_i - u{i+1} = h²f_i
• Siapkan persamaan:

• Yang sama dengan persamaan matriks-vektor:

• Selesaikan persamaan ini dan hasilkan u_i

Salah satu implementasi ini untuk demonstrasi di 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)

Implementasi alternatif tanpa Matriks Aljabar (menggunakan metode 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

Namun Anda dapat menggunakan metode lain untuk menyelesaikan persamaan Laplace. Jika Anda menggunakan metode berulang, Anda harus mengulanginya sampai residual|b-Au|<1e-6 , dengan bmenjadi vektor sisi kananu_L,f_1h²,f_2h²,...

Catatan

Bergantung pada metode solusi Anda, Anda mungkin tidak dapat menyelesaikan contoh dengan tepat untuk solusi yang diberikan. Setidaknya untukN->infinity kesalahan harus mendekati nol.

Celah standar tidak diijinkan , built-in untuk PDE diizinkan.

Bonus

Bonus -30% untuk menampilkan solusi, baik grafis atau ASCII-art.

Kemenangan

Ini codegolf, jadi kode terpendek dalam byte menang!

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

+0,7 byte per komentar @ flawr.

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

Ini memplot output.

misalnya

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

Penjelasan

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

Selesaikan untuk fungsinya u.

Subdivide@#4

Subdivide interval [0,1] menjadi bagian N (input 4).

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

Peta uke output dari Subdivide.

ListPlot[ ... ]

Plot hasil akhirnya.

Solusi non-grafik: 58 byte

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

Matlab, 84, 81.2 79.1 byte = 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)

Perhatikan bahwa dalam contoh ini saya menggunakan vektor baris, ini berarti bahwa matriks Aditransposisikan. fdiambil sebagai pegangan fungsi, a,badalah sisi Dirichlet kiri / kanan yang dikelantang.

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

Sebagai contoh f = 10*x, u_L = 0 u_R = 1, N = 15ini menghasilkan:

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

Sunting: Menyimpan beberapa byte berkat @Angs!

Jika seseorang ingin mengedit gambar, jangan ragu untuk melakukannya. Ini pada dasarnya merupakan adaptasi R dari versi matlab dan python yang diposting.

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

Tidak digabungkan & dijelaskan:

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 
}

Contoh & uji kasus:

Fungsi yang dinamai dan tidak diubah dapat disebut menggunakan:

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)

Perhatikan bahwa fargumennya adalah fungsi-R.

Untuk menjalankan versi golf cukup gunakan (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])

Keterbacaannya cukup terpukul. Tidak Disatukan:

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: Mencetak dalam 83 71 byte.

Biar kulihat:

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

Doh!

Aksioma, 579 460 byte

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 itu dan uji

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)

fungsi untuk pertanyaannya adalah m (,,,) kode di atas dimasukkan ke dalam file "file.input" dan memuat dalam Aksioma. Hasilnya tergantung dari fungsi digit ().

jika seseorang berpikir itu bukan golf => dia dapat menunjukkan bagaimana melakukannya ... terima kasih

PS

tampaknya 6 angka di samping itu. untuk e ^ (x ^ 2) tidak ok di sini atau dalam contoh tetapi di sini saya meningkatkan angka tetapi angka tidak berubah ... bagi saya itu berarti angka dalam contoh salah. Mengapa yang lainnya belum menunjukkan jumlahnya?

untuk dosa (2 *% pi * x) ada masalah juga

"Di sini solusi yang tepat adalah u = (sin (2 * π * x)) / (4 * π ^ 2) +1" "saya menyalin solusi yang tepat untuk x = 1/19:

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

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

1.008224733636964333380661957738992274267070440829381577926...
1.0083001
  1234
1.00822473

1,0083001 diusulkan karena hasilnya berbeda dalam digit ke-4 dari hasil nyata, 1,00822473 ... (dan bukan ke-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