Может ли приближенный якобиан с конечными разностями вызвать нестабильность в методе Ньютона?

13

Я реализовал решатель обратной Эйлера в Python 3 (используя Numpy). Для собственного удобства и в качестве упражнения я также написал небольшую функцию, которая вычисляет аппроксимацию градиента с конечной разностью, чтобы мне не всегда приходилось определять якобиан аналитически (если это вообще возможно!).

Используя описания, приведенные в Ascher и Petzold 1998 , я написал эту функцию, которая определяет градиент в данной точке x:

def jacobian(f,x,d=4):
    '''computes the gradient (Jacobian) at a point for a multivariate function.

    f: function for which the gradient is to be computed
    x: position vector of the point for which the gradient is to be computed
    d: parameter to determine perturbation value eps, where eps = 10^(-d).
        See Ascher und Petzold 1998 p.54'''

    x = x.astype(np.float64,copy=False)
    n = np.size(x)
    t = 1 # Placeholder for the time step
    jac = np.zeros([n,n])
    eps = 10**(-d)
    for j in np.arange(0,n):
        yhat = x.copy()
        ytilde = x.copy()
        yhat[j] = yhat[j]+eps
        ytilde[j] = ytilde[j]-eps
        jac[:,j] = 1/(2*eps)*(f(t,yhat)-f(t,ytilde))
    return jac

Я проверил эту функцию, взяв многомерную функцию для маятника и сравнив символический якобиан с численно определенным градиентом для диапазона точек. Меня порадовали результаты теста, ошибка была около 1е-10. Когда я решил ODE для маятника, используя приближенный якобиан, он также работал очень хорошо; Я не смог обнаружить никакой разницы между ними.

Затем я попытался проверить его с помощью следующего PDE (уравнение Фишера в 1D):

TUзнак равноИкс(КИксU)+λ(U(С-U))

используя конечно-разностную дискретизацию.

Теперь метод Ньютона взрывается с первого шага:

/home/sfbosch/Fisher-Equation.py:40: RuntimeWarning: overflow encountered in multiply
  du = (k/(h**2))*np.dot(K,u) + lmbda*(u*(C-u))
./newton.py:31: RuntimeWarning: invalid value encountered in subtract
  jac[:,j] = 1/(2*eps)*(f(t,yhut)-f(t,yschlange))
Traceback (most recent call last):
  File "/home/sfbosch/Fisher-Equation.py", line 104, in <module>
    fisher1d(ts,dt,h,L,k,C,lmbda)
  File "/home/sfbosch/Fisher-Equation.py", line 64, in fisher1d
    t,xl = euler.implizit(fisherode,ts,u0,dt)
  File "./euler.py", line 47, in implizit
    yi = nt.newton(g,y,maxiter,tol,Jg)
  File "./newton.py", line 54, in newton
    dx = la.solve(A,b)
  File "/usr/lib64/python3.3/site-packages/scipy/linalg/basic.py", line 73, in solve
    a1, b1 = map(np.asarray_chkfinite,(a,b))
  File "/usr/lib64/python3.3/site-packages/numpy/lib/function_base.py", line 613, in asarray_chkfinite
    "array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs

Это происходит для различных значений eps, но, как ни странно, только тогда, когда размер пространственного шага PDE и размер временного шага установлены так, что условие Куранта-Фридрихса-Леви не выполняется. В противном случае это работает. (Это поведение, которое вы ожидаете, если решите с форвардом Эйлером!)

Для полноты вот функция для метода Ньютона:

def newton(f,x0,maxiter=160,tol=1e-4,jac=jacobian):
    '''Newton's Method.

    f: function to be evaluated
    x0: initial value for the iteration
    maxiter: maximum number of iterations (default 160)
    tol: error tolerance (default 1e-4)
    jac: the gradient function (Jacobian) where jac(fun,x)'''

    x = x0
    err = tol + 1
    k = 0
    t = 1 # Placeholder for the time step
    while err > tol and k < maxiter:
        A = jac(f,x)
        b = -f(t,x)
        dx = la.solve(A,b)
        x = x + dx
        k = k + 1
        err = np.linalg.norm(dx)
    if k >= maxiter:
        print("Maxiter reached. Result may be inaccurate.")
        print("k = %d" % k)
    return x

(Функция la.solve является scipy.linalg.solve.)

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

Я вижу в отладчике, что newton () управляет 35 итерациями, прежде чем произойдет ошибка. Это число остается неизменным для каждого EPS, который я пробовал.

Дополнительное наблюдение: когда я вычисляю градиент с помощью FDA и функции, используя начальное условие в качестве входных данных, и сравниваю их при изменении размера эпсилона, ошибка возрастает с уменьшением эпсилона. Я бы ожидал, что сначала он будет большим, затем станет меньше, а затем снова увеличится по мере сжатия эпсилона. Поэтому ошибка в моей реализации якобиана является разумным предположением, но если это так, то она настолько тонка, что я не могу ее увидеть. РЕДАКТИРОВАТЬ: я изменил jacobian () для использования вперед вместо центральных различий, и теперь я наблюдаю ожидаемое развитие ошибки. Однако newton () все еще не сходится. Наблюдая dx в итерации Ньютона, я вижу, что он только растет, нет даже флуктуации: он почти удваивается (коэффициент 1,9) с каждым шагом, причем коэффициент постепенно увеличивается.

Ашер и Петцольд упоминают, что разностные аппроксимации для якобиана не всегда работают хорошо. Может ли приближенный якобиан с конечными разностями вызвать нестабильность в методе Ньютона? Или причина в другом? Как еще я могу подойти к этой проблеме?

Стивен Бош
источник
1
«Я уверен, что моя обратная реализация Эйлера в порядке, потому что я проверил ее, используя функцию для якобиана, и получил стабильные результаты». Просьба уточнить. Вы говорите, что запускаете ту же проблему с точным якобианом, и решение сходится к точному решению PDE? Это важная информация.
Дэвид Кетчон
@DavidKetcheson Да, это то, что я говорю. Извиняюсь, если моя терминология неверна или неполна. (Полагаю, мне следовало бы сказать: «Я получаю стабильные и ожидаемые результаты».)
Стивен Бош

Ответы:

3

Более длинный комментарий, чем все остальное:

Используя описания, приведенные в Ascher и Petzold 1998, я написал эту функцию, которая определяет градиент в данной точке x:

Посмотрите на код разностного приближения SUNDIALS, чтобы получить лучшее представление о том, что вы должны делать в реализации. Ашер и Петцольд - хорошая книга для начала, но SUNDIALS на самом деле используются в производственных работах и, следовательно, были лучше протестированы. (Также SUNDIALS относится к DASPK, над которым работал Петцольд.)

Ашер и Петцольд упоминают, что разностные аппроксимации для якобиана не всегда работают хорошо. Может ли приближенный якобиан с конечными разностями вызвать нестабильность в методе Ньютона?

Эмпирически, приближенные якобианы могут вызывать ошибки сходимости в методе Ньютона. Я не знаю, что я бы охарактеризовал их как «нестабильность»; в некоторых случаях просто невозможно достичь желаемой погрешности в критериях завершения. В других случаях это может проявиться как нестабильность. Я почти уверен, что есть более количественный результат по этому явлению в книге численных методов Хайама или в обсуждении Хайрером и Ваннером W-методов.

Или причина в другом? Как еще я могу подойти к этой проблеме?

Это зависит от того, где вы думаете, ошибка может быть. Если вы очень уверены в своей реализации отсталого Эйлера, я бы не стал там начинать. Опыт сделал меня параноиком в моих реализациях численных методов, поэтому, если бы это был я, я бы начал с кодирования нескольких действительно базовых тестовых задач (пара негибких и жестких линейных задач, уравнения теплопроводности с приближением центрированной конечной разности, и тому подобное), и я бы использовал метод изготовленных решений, чтобы убедиться, что я знаю, каким будет решение, и с чем мне следует сравнивать.

Тем не менее, вы уже сделали кое-что из этого:

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

Это будет следующее, что я бы проверил: используйте аналитический якобиан. После этого вы также можете посмотреть на экстремальные собственные значения вашего якобиана с конечной разницей, если вероятность того, что вы находитесь в нестабильной области отсталого Эйлера, невелика. Изучение экстремальных собственных значений вашего аналитического якобиана в качестве основы для сравнения может дать вам некоторое представление. Предполагая, что все проверены, проблема, вероятно, в решении Ньютона.

Джефф Оксберри
источник
Спасибо за вдумчивый анализ (плюс подсказка SUNDIALS и альтернативные источники). Мой профессор предложил установить lambda = 0, утверждая, что FDA для PDE становится линейной, поэтому мы ожидаем, что FDA Jacobian будет равно аналитическому Jacobian. Когда я делаю это, он управляет тремя временными шагами, ньютон () каждый раз достигает максимума, прежде чем, наконец, взорваться, как раньше.
Стивен Бош
Он также сказал, что не было обычной практики использовать приближенные якобианы для решения PDE, и предположил, что у него могут быть проблемы из-за множества степеней свободы (без объяснения причин, хотя только что взглянул на обсуждение W-методов Хайрером и Ваннером, Я вижу, что это, вероятно, не тривиально).
Стивен Бош
1
Заявление вашего профессора несколько удивительно, учитывая количество литературы по этому вопросу, например, этот знаменитый обзор Кнолля и Кейса . Возможно, мне следовало бы сослаться на эту статью в своем ответе, так как источники в библиографии также могут помочь в диагностике ваших проблем.
Джефф Оксберри