maroba / findiff

Python package for numerical derivatives and partial differential equations in any number of dimensions.
MIT License
420 stars 60 forks source link

funtions with multiple arguments #14

Closed RoyLarson closed 5 years ago

RoyLarson commented 5 years ago

I am trying to use this on a function with multiple arguments like:

from findiff import FinDiff
import numpy as np
def f(x, y):
    return 5*x**2 - 5*x + 10*y**2 -10*y

diff = .01
x = np.linspace(-10,0,11)
y = np.linspace(0,10,11)

d_dx = FinDiff(0, diff)

df_dx = d_dx(f)

This gives me "AttributeError: 'function' object has no attribute 'shape'" Then I tried

df_dx = d_dx(f(x,y))

which produces the result

array([[-1150.],
       [ -850.],
       [ -550.],
       [ -250.],
       [   50.],
       [  350.],
       [  650.],
       [  950.],
       [ 1250.],
       [ 1550.],
       [ 1850.]])

but that is an incorrect result as the values should be from the actual partial derivative

array([[-105.],
       [ -95.],
       [ -85.],
       [ -75.],
       [ -65.],
       [ -55.],
       [ -45.],
       [ -35.],
       [ -25.],
       [ -15.],
       [  -5.]])

Can you help me understand how to use multiple arguments on the function?

Thanks

maroba commented 5 years ago

The problem is that you tried to feed the function with two 1D-arrays of coordinate values, but actually you need a 2D grid in your example. You can create the 2D grid arrays from the coordinate arrays using the numpy.meshgrid function:

from findiff import FinDiff
import numpy as np
def f(x, y):
    return 5*x**2 - 5*x + 10*y**2 -10*y

diff = .01
x = np.linspace(-10,0,11)
y = np.linspace(0,10,11)

X, Y = np.meshgrid(x, y, indexing='ij')

d_dx = FinDiff(0, diff)

df_dx = d_dx(f(X, Y))

np.meshgrid works for any number of dimensions.

FinDiff determines the grid dimension by the arrays that you apply it to. If you pass a 1D-array to a FinDiff object, it thinks that your grid is only 1D. Incidentally, if you had calculated the partial derivative with respect to y instead of x in your original implementation, it would have thrown an IndexError exception, because there is no second index in the 1D array. All in all, if you use grid arrays, all should be fine. :-)

RoyLarson commented 5 years ago

This is the output I expect when I inputted what I put in.

  array([[-105.],
    [ -95.],
    [ -85.],
    [ -75.],
    [ -65.],
    [ -55.],
    [ -45.],
    [ -35.],
    [ -25.],
    [ -15.],
    [  -5.]])

This is the output that I get when I input what you say to

array([[-10500., -10500., -10500., -10500., -10500., -10500., -10500.,
    -10500., -10500., -10500., -10500.],
    [ -9500.,  -9500.,  -9500.,  -9500.,  -9500.,  -9500.,  -9500.,
    -9500.,  -9500.,  -9500.,  -9500.],
    [ -8500.,  -8500.,  -8500.,  -8500.,  -8500.,  -8500.,  -8500.,
    -8500.,  -8500.,  -8500.,  -8500.],
    [ -7500.,  -7500.,  -7500.,  -7500.,  -7500.,  -7500.,  -7500.,
    -7500.,  -7500.,  -7500.,  -7500.],
    [ -6500.,  -6500.,  -6500.,  -6500.,  -6500.,  -6500.,  -6500.,
    -6500.,  -6500.,  -6500.,  -6500.],
    [ -5500.,  -5500.,  -5500.,  -5500.,  -5500.,  -5500.,  -5500.,
    -5500.,  -5500.,  -5500.,  -5500.],
    [ -4500.,  -4500.,  -4500.,  -4500.,  -4500.,  -4500.,  -4500.,
    -4500.,  -4500.,  -4500.,  -4500.],
    [ -3500.,  -3500.,  -3500.,  -3500.,  -3500.,  -3500.,  -3500.,
    -3500.,  -3500.,  -3500.,  -3500.],
    [ -2500.,  -2500.,  -2500.,  -2500.,  -2500.,  -2500.,  -2500.,
    -2500.,  -2500.,  -2500.,  -2500.],
    [ -1500.,  -1500.,  -1500.,  -1500.,  -1500.,  -1500.,  -1500.,
    -1500.,  -1500.,  -1500.,  -1500.],
    [  -500.,   -500.,   -500.,   -500.,   -500.,   -500.,   -500.,
    -500.,   -500.,   -500.,   -500.]])

The numbers are still off, and there are more of them.

Also how do I get just one point

df_0_0_dx = d_dx(f((0,0)))

If you want to look at the jupyter notebook I am trying to figure this out on it is

https://github.com/RoyLarson/Simple_NG_Piping_Pressure_Drop/blob/initial_commit/tests/findiff_test.ipynb

RoyLarson commented 5 years ago

I figured this out. I thought your program took at function evaluated the function along with the given dx and sent back the result at the point. It instead evaluates the finite differences for already evaluated points from a function so the grid had to already be computed.