Численное интегрирование по матрице функций, SymPy и SciPy

Ответ 1

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

А именно, ваш расчет, по-видимому, диагонален в том смысле, что k1 и k2 являются двумя формами k = g^T X g, где X - некоторая матрица 5x5 (с дифференциальными операциями внутри, но это не имеет значения), и g составляет 5xM, при этом M большой. Поэтому k[i,j] = g.T[i,:] * X * g[:,j].

Итак, вы можете просто заменить

for j in xrange(1,n+1):
    for i in xrange(1,m+1):
        g1 += [uu(i,j,x,t),          0,          0,          0,          0]
        g2 += [          0,vv(i,j,x,t),          0,          0,          0]
        g3 += [          0,          0,ww(i,j,x,t),          0,          0]
        g4 += [          0,          0,          0,bx(i,j,x,t),          0]
        g5 += [          0,          0,          0,          0,bt(i,j,x,t)]
g = Matrix( [g1, g2, g3, g4, g5] )

с

i1 = Symbol('i1')
j1 = Symbol('j1')
g1 = [uu(i1,j1,x,t),          0,          0,          0,          0]
g2 = [          0,vv(i1,j1,x,t),          0,          0,          0]
g3 = [          0,          0,ww(i1,j1,x,t),          0,          0]
g4 = [          0,          0,          0,bx(i1,j1,x,t),          0]
g5 = [          0,          0,          0,          0,bt(i1,j1,x,t)]
g_right = Matrix( [g1, g2, g3, g4, g5] )

i2 = Symbol('i2')
j2 = Symbol('j2')
g1 = [uu(i2,j2,x,t),          0,          0,          0,          0]
g2 = [          0,vv(i2,j2,x,t),          0,          0,          0]
g3 = [          0,          0,ww(i2,j2,x,t),          0,          0]
g4 = [          0,          0,          0,bx(i2,j2,x,t),          0]
g5 = [          0,          0,          0,          0,bt(i2,j2,x,t)]
g_left = Matrix( [g1, g2, g3, g4, g5] )

и

tmp = evaluateExpr( B*g )
k1 = r*tmp.transpose() * F * tmp
k2 = r*g.transpose()*evaluateExpr(Bc*g)
k2 = evaluateExpr( k2 )

по

tmp_right = evaluateExpr( B*g_right )
tmp_left = evaluateExpr( B*g_left )
k1 = r*tmp_left.transpose() * F * tmp_right
k2 = r*g_left.transpose()*evaluateExpr(Bc*g_right)
k2 = evaluateExpr( k2 )

Не тестировал (прошлый час), но вы поняли.

Теперь вместо того, чтобы иметь огромную символическую матрицу, которая делает все медленным, у вас есть два матричных индекса для индексов пробной функции и свободные параметры i1,j1 и i2,j2, которые играют свою роль, и вы должны подставить в них целые числа в конец.

Так как матрица для lambdify равна только 5x5 и должна быть окутана только один раз за пределами всех циклов, накладные расходы на лямблизацию и упрощение исчезнут. Более того, проблема легко помещается в память даже при больших т, п.

Интеграция не выполняется быстрее, но поскольку выражения очень малы, вы можете легко, например, сбрасывать их в Fortran или делать что-то еще умное.

Ответ 2

quadpy (мой проект) делает векторизованное числовое интегрирование. Это

from numpy import sin, cos, pi
import quadpy


def f(X):
    x, t = X
    return [
        [(sin(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*cos(t)*x)*cos(3-0.1*x)*cos(t)],
        [(cos(2-0.1*x)*sin(t)*x+sin(2-0.1*x)*cos(t)*x)*sin(3-0.1*x)*cos(t)],
        [(cos(2-0.1*x)*sin(t)*x+cos(2-0.1*x)*sin(t)*x)*sin(3-0.1*x)*sin(t)]
        ]


x1 = 30
x2 = 75
t1 = 0
t2 = 2*pi

sol = quadpy.quadrilateral.integrate(
        f,
        [[x1, t1], [x2, t1], [x2, t2], [x1, t2]],
        quadpy.quadrilateral.Product(quadpy.line_segment.GaussLegendre(5))
        )

print(sol)

дает

[[ 1456.3701526 ]
 [ 2620.60490653]
 [ 5034.5831071 ]]

Тайминги:

%timeit quadpy.quadrilateral.integrate(f, [[x1, t1], [x2, t1], [x2, t2], [x1, t2]], q)
1000 loops, best of 3: 219 µs per loop

Это приводит к резкому ускорению в загружаемом примере:

import numpy
array2mat = [{'ImmutableMatrix': numpy.array}, 'numpy']
k1_lambda = lambdify( (x,t), k1, modules=array2mat)
print 'Finished lambdifying k1:', time.clock()
import quadpy
sol = quadpy.quadrilateral.integrate(
    lambda X: k1_lambda(X[0], X[1]),
    [[x1, t1], [x2, t1], [x2, t2], [x1, t2]],
    quadpy.quadrilateral.Product(quadpy.line_segment.GaussLegendre(5))
    )

Вывод:

Start: 0.040001
Finished trial functions: 0.379929
Finished evaluating differential equations: 2.669536
Finished lambdifying k1: 29.961808
Finished integrating k1: 30.106988
Finished lambdifying and integrating k2: 34.229007
Finished calculating eigenvalues and eigenvectors: 34.229924

Обратите внимание, что quadpy не выполняет адаптивную квадратуру, поэтому разумно выбирайте вашу схему.