Метод стрельбы в единичных матричных предупреждениях - PullRequest
Новая
       42

Метод стрельбы в единичных матричных предупреждениях

0 голосов
/ 04 марта 2019

Я пытаюсь реализовать метод стрельбы в python и получаю сообщения об ошибках для единичных матриц.Я хочу решить xy''' = -y'y+x^2-1 с граничными условиями y(1)=1, y(2)=1, y'(2)=4.Я преобразовал это в уравнения первого порядка M (x, y ) y '= f (x, y ), с: 

M = [1  0  0]
    [0  1  0]
    [0  0  x]

**y** = [ y ]
        [ y']
        [y''] 
f = [     y'   ]
    [     y''  ]
    [-y'y+x^2-1]

Я реализовал метод стрельбы, как описано в Числовых рецептах 18.1.Общая процедура:

  1. Укажите граничные условия как с левой, так и с правой стороны.В моем случае это y (1) = 1, y (2) = 1, y '(2) = 4.Укажите начальное предположение для вектора (y, y ', y' ').Я выбрал (1, 1, 2).Первый элемент фиксирован, так как y (1) = 1, но два других элемента составляют свободно определяемый вектор V , содержащий другие граничные условия на LHS.
  2. Интегрируем уравненияот х = 1 до х = 2.Определите вектор несоответствия F , рассчитанный как разность между предписанными граничными условиями (y (2) = 1 и y '(2) = 4) и значениями, вычисленными посредством интегрирования.
  3. Используйте Ньютона-Рафсона, чтобы найти V , нули F .Выполните это, решив J .δ V = -F и добавив поправку обратно к V . J - это якобиан, найденный численно на каждом временном шаге.

Если я предполагаю, что M обратимо, то это относительно просто решить и работает как ожидалось.Тем не менее, я столкнулся с двумя проблемами и был бы признателен за помощь в их решении.

  1. Если я установлю xrange от -1 до 1, я получу предупреждение о единичной матрице.Я думаю, это потому, что я предположил, что М обратим, но это не так, когда х = 0.Такая же проблема возникает часто, если M имеет какую-либо зависимость от y Это похоже на дифференциальное алгебраическое уравнение.Есть ли какие-нибудь удобные способы решить эту проблему в python?Не похоже, что есть какие-то встроенные в scipy.
  2. Если я изменю свое граничное условие при x = 2 на y (2) = - 1, я также получу предупреждение о единичной матрице.Я не уверен, почему это происходит.Может быть, это просто неспособность дифференциальных уравнений удовлетворять граничным условиям или что-то не так с моим кодом?Если первое, есть ли способ предсказать это заранее?

Предыдущий ответ предположил, что лучшее предположение о начальных условиях может помочь во втором единственном случае, ноЯ не уверен, как улучшить мои предположения без большого количества грубой силы.

Вот мой код: 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

### Solving the differential equation xy''' = -y'y+x^2-1  
### if x range doesn't include 0, then M is invertible, so equations become:
def fn(y,x):
    return np.array([y[1],y[2],(-y[0]*y[1]+x**2-1)/x])

### Boundary values of y(1)=1, y(2)=1, y'(2)=4

x = np.linspace(1.,2.,101)
y0 = np.array([1., 1., 2.]) ### Only the first of these is fixed. The others are variable
sol = np.array([1.,4.]) ### Specify y(2)=1 y'(2)=4
Delta=1
max_itrs = 100

eps = 1e-3

def shoot(y0, Delta, fn, x, sol, eps, max_itrs):
    ### Throughtout this function, the indexing only works with our specific boundary conditions
    ### But can be easily modified or generalized to accommodate other cases
    for i in range(max_itrs):
        F = odeint(fn, y0, x)[-1][:2]-sol ### Define the F "discrepancy" vector at x = 2
        if np.linalg.norm(F)>eps:
            y0_yp_mod = y0.copy()
            y0_yp_mod[1] += Delta ### Augment y'(1) by Delta
            F_yp_mod = odeint(fn, y0_yp_mod, x)[-1][:2]-sol 

            y0_ypp_mod = y0.copy()
            y0_ypp_mod[2] += Delta ### Augment y''(1) by Delta
            F_ypp_mod = odeint(fn, y0_ypp_mod, x)[-1][:2]-sol

            ### Compute the Jacobian
            j11 = (F_yp_mod[0]-F[0])/Delta
            j12 = (F_yp_mod[1]-F[1])/Delta
            j21 = (F_ypp_mod[0]-F[0])/Delta
            j22 = (F_ypp_mod[1]-F[1])/Delta
            jac = np.array([[j11,j12],[j21,j22]])

            ### Solve for delta V:
            d_V = np.linalg.solve(jac,-F)

            ### Augment our initial conditions by delta V
            y0 = np.array([y0[0],y0[1]+d_V[0],y0[2]+d_V[1]])
        else:
            return y0
    print("Failed to achieve convergence after "+str(max_itrs)+" iterations.")
    return y0

y0 = shoot(y0, Delta, fn, x, sol, eps, max_itrs)
Y = odeint(fn, y0, x)
plt.plot(x,Y[:,0],c="r",label="y")
plt.plot(x,Y[:,1],c="g",label="y'")
plt.plot(x,Y[:,2],c="b",label="y''")
plt.scatter([x[0]],y0[0],marker="x",c="r")
plt.scatter([x[-1]],sol[0],marker="x",c="r")
plt.scatter([x[-1]],sol[1],marker="x",c="g")
plt.legend()

Здесь - это выход, который ведет себякак и ожидалось output

...