Résoudre l'équation de Laplace

13

Introduction aux mathématiques numériques

Ceci est le "Bonjour, monde!" des PDE (Partial Differential Equations). L'équation de Laplace ou de diffusion apparaît souvent en physique, par exemple l'équation de la chaleur, la déformation, la dynamique des fluides, etc. et pas chanter "99 bouteilles de bière, ..." cette tâche est donnée en 1D. Vous pouvez interpréter cela comme une robe en caoutchouc attachée à un mur aux deux extrémités avec une certaine force appliquée.

Sur un [0,1]domaine, trouvez une fonction upour une fonction source donnée fet des valeurs limites u_Let u_Rtelles que:

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

u'' désigne la dérivée seconde de u

Cela peut être résolu purement théorique mais votre tâche est de le résoudre numériquement sur un domaine discrétisé x pour les Npoints:

  • x = {i/(N-1) | i=0..N-1}ou 1:{(i-1)/(N-1) | i=1..N}
  • h = 1/(N-1) est l'espacement

Contribution

  • f comme fonction ou expression ou chaîne
  • u_L, en u_Rtant que valeurs à virgule flottante
  • N comme entier> = 2

Production

  • Array, List, une sorte de chaîne séparée de utelle sorte queu_i == u(x_i)

Exemples

Exemple 1

Entrée: f = -2, u_L = u_R = 0, N = 10(Ne pas prendre f=-2mal, ce n'est pas une valeur , mais une fonction constante que le rendement -2pour tous xIl est comme une force de gravité constante sur notre corde..)

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

Il existe une solution exacte simple: u = -x*(1-x)

Exemple 2

Entrée: f = 10*x, u_L = 0 u_R = 1, N = 15(Ici il y a beaucoup de vent debout sur le côté droit)

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

La solution exacte pour cela indique: u = 1/3*(8*x-5*x^3)

Exemple 3

Entrée: f = sin(2*pi*x), u_L = u_R = 1, N = 20(Quelqu'un a brisé la gravité ou il y a une sorte de haut et sous le vent)

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

Ici, la solution exacte est u = (sin(2*π*x))/(4*π^2)+1

Exemple 4

Entrée: f = exp(x^2), u_L = u_R = 0,N=30

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

Notez la légère asymétrie

FDM

Une méthode possible pour résoudre ce problème est la méthode des différences finies :

  • réécrire -u_i'' = f_icomme
  • (-u_{i-1} + 2u_i - u{i+1})/h² = f_i ce qui équivaut
  • -u_{i-1} + 2u_i - u{i+1} = h²f_i
  • Configurez les équations:

  • Qui sont égaux à une équation matrice-vecteur:

  • Résolvez cette équation et sortez le u_i

Une implémentation de ceci pour démonstration en 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)

Implémentation alternative sans algèbre matricielle (en utilisant la méthode 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

Vous pouvez cependant utiliser toute autre méthode pour résoudre l'équation de Laplace. Si vous utilisez une méthode itérative, vous devez itérer jusqu'à ce que le résidu |b-Au|<1e-6, avec bêtre le vecteur de droiteu_L,f_1h²,f_2h²,...

Remarques

En fonction de votre méthode de résolution, il se peut que vous ne résolviez pas les exemples exactement aux solutions données. Au moins pour N->infinityl'erreur devrait approcher de zéro.

Les failles standard sont interdites , les fonctions intégrées pour les PDE sont autorisées.

Prime

Un bonus de -30% pour l'affichage de la solution, graphique ou ASCII.

Gagnant

Ceci est codegolf, donc le code le plus court en octets gagne!

code-golf math number Karl Napf
la source
Je recommande d'ajouter un exemple qui n'est pas résolu analytiquement, par exemple avec f(x) = exp(x^2).
flawr
@flawr Bien sûr, il a une solution mais la fonction d'erreur est impliquée.
Karl Napf
1
Désolé, ce n'était peut-être pas la bonne expression, est-ce que "l'antidérivative non élémentaire" pourrait être mieux adaptée? Je veux dire des fonctions comme log(log(x))ou sqrt(1-x^4)qui ont une intégrale, qui n'est cependant pas exprimable dans les fonctions élémentaires.
flawr
@flawr Non ça va, la fonction d'erreur n'est pas élémentaire, je voulais juste dire qu'il existe un moyen d'exprimer la solution analytiquement mais u(x) = 1/2 (-sqrt(π) x erfi(x)+sqrt(π) erfi(1) x+e^(x^2)-e x+x-1)n'est pas exactement calculable.
Karl Napf
Pourquoi répéter jusqu'à 1e-6 et ne pas répéter jusqu'à 1e-30?
RosLuP

Réponses:

4

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

+0,7 octets par commentaire de @flawr.

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

Cela trace la sortie.

par exemple

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

entrez la description de l'image ici

Explication

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

Résolvez la fonction u.

Subdivide@#4

Subdivide l'intervalle [0,1] en N (4e entrée) parties.

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

Mappez uà la sortie de Subdivide.

ListPlot[ ... ]

Tracez le résultat final.

Solution non graphique: 58 octets

u/@Subdivide@#4/.NDSolve[-u''@x==#&&u@0==#2&&u@1==#3,u,x]&
JungHwan Min
la source
Cela ne fonctionne pas pourf(x) = exp(x^2)
flawr
Vous voudrez peut-être utiliser NDSolvepour le cas général des solutions non élémentaires.
flawr
6

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

Notez que dans cet exemple, j'utilise des vecteurs de ligne, cela signifie que la matrice Aest transposée. fest considéré comme une poignée de fonction, a,bsont les contraintes de Dirichlet gauche / droite.

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

Pour l'exemple, f = 10*x, u_L = 0 u_R = 1, N = 15cela se traduit par:

flawr
la source
3

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

Edit: enregistré une poignée d'octets grâce à @Angs!

Si quelqu'un souhaite modifier les images, n'hésitez pas à le faire. Il s'agit essentiellement d'une adaptation R des versions matlab et python publiées.

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

Non golfé et expliqué:

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 
}

Exemples et cas de test:

La fonction nommée et non golfée peut être appelée en utilisant:

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)

Notez que l' fargument est une fonction R.

Pour exécuter la version golfée, utilisez simplement (function(...){...})(args)

entrez la description de l'image ici entrez la description de l'image ici entrez la description de l'image ici entrez la description de l'image ici

Billywob
la source
Je pense que vous pouvez laisser tomber le is.numeric(f)chèque si vous déclarez fcomme fonction, il n'est pas nécessaire de le passer directement dans l'appel de fonction au solveur.
Karl Napf
Ah je vois, je ne savais pas qu'il y avait une différence entre ces deux. Eh bien, s'il est plus court, vous pouvez modifier votre solveur pour l'accepter fcomme une fonction afin que vous n'ayez pas à vérifier si c'est une constante (fonction).
Karl Napf
1
@Billywob, il n'est pas nécessaire fd'être toujours numérique. f = (function(x)-2)fonctionne pour le premier exemple, donc il n'y a jamais besoin de le faire rep.
Angs
Vous pouvez utiliser x<-0:10/10;f<-function(x){-2};10^-2*sapply(x,f)si f (x) n'est pas garanti pour être vectorisé ou tout simplement 10^-2*f(x)s'il fest vectorisé ( laplace(Vectorize(function(x)2),0,0,10)
Angs
N'utilisez pas eval, donnez fcomme fonction appropriée.
Angs
2

Haskell, 195168 octets

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

La lisibilité a pris un coup. Non golfé:

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: Impression en 83 71 octets.

Laisse-moi voir:

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

Oh!

Angs
la source
Je ne sais pas grand chose sur Haskell, mais peut-être que la solution sans algèbre matricielle pourrait être plus courte, j'ai ajouté un deuxième exemple d'implémentation.
Karl Napf
@KarlNapf ne se rapproche pas très près Ce n'est que semi-golfé mais il doit utiliser beaucoup de fonctions intégrées verbeuses. Avec l'algèbre matricielle, la majeure partie du code consiste à créer la matrice (64 octets) et l'importation (29 octets). Le résidu et le jacobi prennent beaucoup de place.
Angs
Eh bien, tant pis mais ça valait le coup d'essayer.
Karl Napf
1

Axiom, 579 460 octets

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)

dé-golfer et tester

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)

la fonction pour la question est m (,,,) le code ci-dessus est mis dans le fichier "file.input" et chargé dans Axiom. Le résultat dépend de la fonction digits ().

si quelqu'un pense que ce n'est pas du golf => il ou elle peut montrer comment le faire ... merci

PS

il semble que les 6 chiffres après le. pour e ^ (x ^ 2) ne sont pas corrects ici ou dans les exemples mais ici j'augmente les chiffres mais les nombres ne changent pas ... pour moi cela signifie que les nombres dans l'exemple sont faux. Pourquoi tous les autres n'ont pas montré leurs chiffres?

pour le péché (2 *% pi * x) il y a aussi des problèmes

"Ici, la solution exacte est u = (sin (2 * π * x)) / (4 * π ^ 2) +1" j'ai copié la solution exacte pour x = 1/19: 

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

dans WolframAlpha https://www.wolframalpha.com/input/?i=(sin(2% CF% 80% 2F19))% 2F (4 % CF% 80% 5E2)% 2B1 il en résulte

1.008224733636964333380661957738992274267070440829381577926...
1.0083001
  1234
1.00822473

1.0083001 proposé comme résultat est différent au 4ème chiffre du résultat réel 1.00822473 ... (et non pas au 6ème) 

-- 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
la source
La solution numérique diffère de la solution exacte car le FDM n'est ici que de second ordre, ce qui signifie que seuls les polynômes jusqu'à l'ordre 2 peuvent être représentés exactement. Donc, seul l' f=-2exemple a une solution analytique et numérique correspondante.
Karl Napf
ici, la solution numérique semble correcte si je change les chiffres en 80 ou 70 -> g (sin (2 *% pi * x), 1,1,1 / 19) 1,0082247336 3696433338 0661957738 9922742670 7044082938 1577926908 950765832 l'autre numéro 1.0082247336 3696433338 0661957038 9922 7044082938 1577926 ...
RosLuP