Подгонка функции f (x, y, z) с квадратичным c полиномом

Question

Я пытаюсь согласовать функцию f (x, y, z) со следующим квадратиком c полином:

3-й полином

Некоторые искажены сферическая поверхность в трех измерениях. Проблема связана с вычислением эффективных масс в solid физике состояний.

Вот изображение данных, показывающее, что оно действительно параболически падает во всех направлениях, даже если кривизна в z- направление довольно низкое: 3d параболы

Меня интересуют коэффициенты, которые соответствуют эффективным массам. У меня есть массив координат XYZ, который является регулярным и центрирован на максимуме:

[[ 0.          0.          0.        ]
 [ 0.          0.          0.01282017]
 [ 0.          0.          0.02564034]
 ...
 [-0.05026321 -0.05026321 -0.03846052]
 [-0.05026321 -0.05026321 -0.02564034]
 [-0.05026321 -0.05026321 -0.01282017]]

и соответствующий одномерный массив скалярных значений, по одному для каждой точки. Число точек данных вокруг этого максимума может варьироваться от 100 до 1000.

Это код, который я сейчас пытаюсь использовать для подгонки:

def func(data, mxx, mxy, mxz, myy, myz, mzz):
    x = data[:, 0]
    y = data[:, 1]
    z = data[:, 2]
    return (
        (1 / (2 * mxx)) * (x ** 2)
        + (1 / (1 * mxy)) * (x * y)
        + (1 / (1 * mxz)) * (x * z)
        + (1 / (2 * myy)) * (y ** 2)
        + (1 / (1 * myz)) * (y * z)
        + (1 / (2 * mzz)) * (z ** 2)
    ) + f(0, 0, 0)

energy = data[:, 3]
guess = (mxx, mxy, mxz, myy, myz, mzz)
params, pcov = scipy.optimize.curve_fit(
    func, data, energy, p0=guess, method="trf"
)

Где f (0,0 , 0) это значение функции в (0, 0, 0), которое я получаю с помощью функции scipy.interpolate.griddata.

Для этой задачи массы должны быть отрицательными и иметь значения между - 0,2 и -2, грубо говоря. Я создаю предположительные значения с помощью конечно-разностного дифференцирования.

Однако я не получаю никаких значимых результатов от scipy.interpolate.curve_fit - обычно коэффициенты заканчиваются огромными числами (например, 1e9). Я полностью потерян в этот момент.

Что я делаю не так :(?

Answer 1

Одна из проблем в том, что вы подходите 1/m. Хотя это правильно с точки зрения физики, это плохо с точки зрения алгоритма. Если алгоритму подбора необходимо изменить знак для значений m около нуля, коэффициенты расходятся. Следовательно, лучше подгонять mI = 1/m и делать соответствующие ошибки позже. Здесь я использую leastsq, что требует дополнительных вычислений для ковариационной матрицы (поскольку она возвращает уменьшенную форму). Я подгоняю g() и инверсные массы, но вы можете немедленно воспроизвести свои проблемы, введя f() и непосредственно подгоняя m s.

Второй момент заключается в том, что данные имеют смещение Т.е. если x = y = z = 0, то данные v= -0.0195 Это необходимо ввести в модель.

Наконец, я бы сказал, что у вас уже есть не параболическое поведение c в ваших данных.

Тем не менее, вот как это выглядит:

import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(linewidth=300)

from scipy.optimize import leastsq
from scipy.optimize import curve_fit

data = np.loadtxt( "silicon.csv", delimiter=',' )

def f( x, y, z, mxx, mxy, mxz, myy, myz, mzz, offI ):
    out = 1./(2 * mxx) * x * x
    out += 1./( mxy ) * x * y
    out += 1./( mxz ) * x * z
    out += 1./( 2 * myy ) * y * y
    out += 1./( myz ) * y * z
    out += 1./( 2 * mzz ) * z * z
    out += 1./offI
    return out

def g( x, y, z, mxxI, mxyI, mxzI, myyI, myzI, mzzI, off ):
    out = mxxI / 2 * x * x
    out += mxyI  * x * y
    out += mxzI * x * z
    out += myyI / 2 * y * y
    out += myzI  * y * z
    out += mzzI / 2  * z * z
    out += off
    return out


def residuals( params, indata ):
    out = list()
    for x, y, z, v in indata:
        out.append( v - g( x,y, z, *params ) )
    return out

sol, cov, info, msg, ier = leastsq( residuals,  7*[0], args=( data, ), full_output=True)

s_sq = sum( [x**2 for x in residuals( sol, data) ] )/ (len( data ) - len( sol ) )
print "solution"
print sol

masses = [1/x for x in sol]
print "masses:"
print masses

print "covariance matrix:"
covMX = cov * s_sq
print covMX

print "sum of residuals"
print sum( residuals( sol, data) )

### plotting the cuts

fig = plt.figure('cuts')
ax = dict()
for i in range( 1, 10 ):
    ax[i] = fig.add_subplot( 3, 3, i )

dl = np.linspace( -.2, .2, 25)

#### xx
xdata = [ [ x, v ] for x,y,z,v in data if ( abs(y)<1e-3 and abs(z) < 1e-3 ) ]
vl = np.fromiter( ( f( x, 0, 0, *masses ) for x in dl ), np.float )
ax[1].plot( *zip(*sorted( xdata ) ), ls='', marker='o')
ax[1].plot( dl, vl )

#### xy
xydata = [ [ x, v ] for x, y, z, v in data if ( abs( x - y )<1e-2 and abs(z) < 1e-3 ) ]
vl = np.fromiter( ( f( xy, xy, 0, *masses ) for xy in dl ), np.float )
ax[2].plot( *zip(*sorted( xydata ) ), ls='', marker='o')
ax[2].plot( dl, vl )

#### xz
xzdata = [ [ x, v ] for x, y, z, v in data if ( abs( x - z )<1e-2 and abs(y) < 1e-3 ) ]
vl = np.fromiter( ( f( xz, 0, xz, *masses ) for xz in dl ), np.float )
ax[3].plot( *zip(*sorted( xzdata ) ), ls='', marker='o')
ax[3].plot( dl, vl )

#### yy
ydata = [ [ y, v ] for x, y, z, v in data if ( abs(x)<1e-3 and abs(z) < 1e-3 ) ]
vl = np.fromiter( ( f( 0, y, 0, *masses ) for y in dl ), np.float )
ax[5].plot( *zip(*sorted( ydata ) ), ls='', marker='o' )
ax[5].plot( dl, vl )

#### yz
yzdata = [ [ y, v ] for x, y, z, v in data if ( abs( y - z )<1e-2 and abs(x) < 1e-3 ) ]
vl = np.fromiter( ( f( 0, yz, yz, *masses ) for yz in dl ), np.float )
ax[6].plot( *zip(*sorted( yzdata ) ), ls='', marker='o')
ax[6].plot( dl, vl )

#### zz
zdata = [ [ z, v ] for x, y, z, v in data if ( abs(x)<1e-3 and abs(y) < 1e-3 ) ]
vl = np.fromiter( ( f( 0, 0, z, *masses ) for z in dl ), np.float )
ax[9].plot( *zip(*sorted( zdata ) ), ls='', marker='o' )
ax[9].plot( dl, vl )

#### some diag
ddata = [ [ z, v ] for x, y, z, v in data if ( abs(x - y)<1e-3 and abs(x - z) < 1e-3 ) ]
vl = np.fromiter( ( f( d, d, d, *masses ) for d in dl ), np.float )
ax[7].plot( *zip(*sorted( ddata ) ), ls='', marker='o' )
ax[7].plot( dl, vl )

#### some other diag
ddata = [ [ z, v ] for x, y, z, v in data if ( abs(x - y)<1e-3 and abs(x + z) < 1e-3 ) ]
vl = np.fromiter( ( f( d, d, -d, *masses ) for d in dl ), np.float )
ax[8].plot( *zip(*sorted( ddata ) ), ls='', marker='o' )
ax[8].plot( dl, vl )

plt.show()

Это дает следующий вывод:

solution
[-1.46528595  0.25090717  0.25090717 -1.46528595  0.25090717 -1.46528595 -0.01993436]
masses:
[-0.6824606499739905, 3.985537743156507, 3.9855376943660676, -0.6824606473928339, 3.9855377322848344, -0.6824606467055248, -50.16463861555409]
covariance matrix:
[
[ 4.76417852e-03 -1.46907683e-12 -8.57639600e-12 -2.21281938e-12 -2.38444957e-12  8.42981521e-12 -2.70034183e-05]
[-1.46907683e-12  9.17104397e-04 -7.10573582e-13  1.32125214e-11  7.44553140e-12  1.29909935e-11 -1.11259046e-13]
[-8.57639600e-12 -7.10573582e-13  9.17104389e-04 -8.60004172e-12 -6.14797647e-12  8.27070243e-12  3.11127064e-14]
[-2.21281914e-12  1.32125214e-11 -8.60004172e-12  4.76417860e-03 -4.20477032e-12  9.20893224e-12 -2.70034186e-05]
[-2.38444957e-12  7.44553140e-12 -6.14797647e-12 -4.20477032e-12  9.17104395e-04  1.50963408e-11 -7.28889534e-14]
[ 8.42981530e-12  1.29909935e-11  8.27070243e-12  9.20893175e-12  1.50963408e-11  4.76417849e-03 -2.70034182e-05]
[-2.70034183e-05 -1.11259046e-13  3.11127064e-14 -2.70034186e-05 -7.28889534e-14 -2.70034182e-05  5.77019926e-07]
]
sum of residuals
4.352727352163743e-09

... и здесь некоторые 1d срезы, которые показывают некоторое существенное отклонение от поведения параболи c, если один не находится на одной из главных осей.