Failed to interpolate Gaussian distribution

1

Hi all,

I have a question regarding interpolation. I wanted to interpolate a Gaussian distribution for a boundary condition. However, I could not obtain a correct interpolation of the expression defining the distribution. I came across a similar issue when interpolating a function defined by exponential function. Do you have any ideas on what is going on here? Note that the version of Dolphin is 2019.2.0.

import dolfin as dl
import numpy as np
import matplotlib.pyplot as plt
N = 100 # number of elements
T = 10
mesh = dl.IntervalMesh(N, 0, T)
degree = 1
V = dl.FunctionSpace(mesh, ā€œLagrangeā€, degree)
mu = 5
sigma = 1.0
u_ini = dl.Expression(ā€˜1/(sigma* sqrt(2* pi))* exp(-1/2* (pow((x[0] - mu)/sigma, 2)))ā€™, mu = mu, sigma = sigma, degree = 1)
u_n = dl.interpolate(u_ini, V)

time = np.linspace(0, T, 101)

u_true = 1/(sigma* np.sqrt(2* 3.14))* np.exp(-1/2* ((time - mu)/sigma)**2)

dl.plot(u_n)
plt.plot(time, u_true)
plt.ylim(-0.1, 1.1)

I appreciate your help.

2

The problem is in your definition of u_ini. You should be careful of order of evaluation. Whichever of * or / lies on the left of the other gets evaluated first, and always before + and -. For this simple case, you should be able to see by evaluating the below expression and check that

(pow((x[0] - mu)/sigma, 2.0))*1/2

evaluates to your desired result, whereas

1/2*(pow((x[0] - mu)/sigma, 2.0))

evaluates to zero (0). And exp(0) always evaluates to 1 no matter what your x[0] is. You can explicitly specify floating point division as a default in your expression (1./2 vs 1/2) by

dl.Expression('1/(sigma*sqrt(2*pi))*exp(-1./2*(pow((x[0] - mu)/sigma, 2.0)))', mu = mu, sigma = sigma, degree = 1)

Also, you can enclose your code in triple backticks ``` to ensure proper readability and execution.

3 Likes
3

Hi Bhavesh,

Thank you for your help! I did not notice the problem is in the floating division since I am not familiar with C++ syntax (so used to Python easy oneā€¦). Thank you so much!

1 Like
4

thanks for the awesome information.

5

thanks my issue has been fixed.

6

Hello @dokken

How can I make a graph of the distribution calculated through the expression?

7

Project the expression into an appropriate function space, and then plot it with either matplotlib or paraview

8

This is my code. My mesh is in 2D vertical.

from dolfin import *
import numpy as np

L = 100
mesh = UnitSquareMesh(L, L)
    
V = FunctionSpace(mesh, "Lagrange", 1)
    
L = 11000
mu_phi = 0.5
x0_phi = 0.5 * L
phi0 = Expression('exp(-mu_phi*pow(x[0] - x0_phi, 2.0))', degree=1, mu_phi=mu_phi, x0_phi=x0_phi)
phia = interpolate(phi0, V)

file = File("phi/phi0.pvd")
file << phia
9

I am not sure what you expect. Just plotting this with wolfram alpha gives you 0 everywhere. Due to your large L (are you sure L shouldnā€™t be 1?): plot e^(-0.5*(x-5500)**2) for x in 0, 1 - Wolfram|Alpha

10

Sorry, I had a 11000 long mesh but for the example I used a fenics mesh. I just want to know if I did the distribution right by looking at the graph. But I canā€™t do it with the fenics plot nor with the paraview.

from dolfin import *
import numpy as np

mesh = RectangleMesh(Point(0.0, 0.0), Point(100.0, 50.0), 10, 10)

V = FunctionSpace(mesh, "Lagrange", 1)
mu_phi = 0.5
L = 100

V = FunctionSpace(mesh, "Lagrange", 1)

mu_phi = 0.5
x0_phi = 0.5 * L
phi0 = Expression('exp(-mu_phi*pow(x[0] - x0_phi, 2.0))', degree=1, mu_phi=mu_phi, x0_phi=x0_phi)
phia = interpolate(phi0, V)

file = File("phi/phi0.pvd")
file << phia