Конечный элементный код начинающего не решает уравнение должным образом

Question

Я пытаюсь написать код для решения чрезвычайно сложного дифференциального уравнения: x' = 1 с помощью метода конечных элементов. Насколько я понял, я могу получить решение u как

с базисными функциями phi_i(x), в то время как я могу получить u_i как решение системы линейных уравнений:

с дифференциальным оператором D (здесь только первая производная). В качестве основы я использую функцию палатки:

    def tent(l, r, x):
    m = (l + r) / 2 
    if x >= l and x <= m:
        return (x - l) / (m - l)
    elif x < r and x > m:
        return (r - x) / (r - m)
    else:
        return 0

def tent_half_down(l,r,x):
    if x >= l and x <= r:
        return (r - x) / (r - l)
    else:
        return 0

def tent_half_up(l,r,x):
    if x >= l and x <= r:
        return (x - l) / (r - l)
    else:
        return 0

def tent_prime(l, r, x):
    m = (l + r) / 2 
    if x >= l and x <= m:
        return 1 / (m - l)
    elif x < r and x > m:
        return 1 / (m - r)
    else:
        return 0

def tent_half_prime_down(l,r,x):
    if x >= l and x <= r:
        return - 1 / (r - l)
    else:
        return 0

def tent_half_prime_up(l, r, x):
    if x >= l and x <= r:
        return 1 / (r - l)
    else:
        return 0

def sources(x):
    return 1

Дискретизация моего пространства:

n_vertex = 30
n_points = (n_vertex-1) * 40
space = (0,5)
x_space = np.linspace(space[0],space[1],n_points)
vertx_list = np.linspace(space[0],space[1], n_vertex)
tent_list = np.zeros((n_vertex, n_points))
tent_prime_list = np.zeros((n_vertex, n_points))

tent_list[0,:] = [tent_half_down(vertx_list[0],vertx_list[1],x) for x in x_space]
tent_list[-1,:] = [tent_half_up(vertx_list[-2],vertx_list[-1],x) for x in x_space]
tent_prime_list[0,:] = [tent_half_prime_down(vertx_list[0],vertx_list[1],x) for x in x_space]
tent_prime_list[-1,:] = [tent_half_prime_up(vertx_list[-2],vertx_list[-1],x) for x in x_space]
for i in range(1,n_vertex-1):
    tent_list[i, :] = [tent(vertx_list[i-1],vertx_list[i+1],x) for x in x_space]
    tent_prime_list[i, :] = [tent_prime(vertx_list[i-1],vertx_list[i+1],x) for x in x_space]

Расчет системы линейных уравнений:

b = np.zeros((n_vertex))
A = np.zeros((n_vertex,n_vertex))
for i in range(n_vertex):
    b[i] = np.trapz(tent_list[i,:]*sources(x_space))
    for j in range(n_vertex):
        A[j, i] = np.trapz(tent_prime_list[j] * tent_list[i])

А затем решить и восстановить его

u = np.linalg.solve(A,b)
sol = tent_list.T.dot(u)

Но это не работает, я только получаю некоторые взлеты и падения. Что я делаю неправильно?

Answer 1

Сначала пара комментариев по терминологии и обозначениям:

1) Вы используете слабую формулировку, хотя вы сделали это неявно. Слабая формулировка не имеет ничего общего с порядком участвующих производных. Он слабый, потому что вы не удовлетворяете дифференциальному уравнению точно в каждом месте. FE минимизирует взвешенный остаток решения, интегрированного в домен. Функции phi_j фактически дискретизируют весовую функцию. Различие, когда у вас есть только производные первого порядка, состоит в том, что вам не нужно применять теорему о расходимости Гаусса (которая упрощает интеграцию по частям для одного измерения), чтобы исключить производные второго порядка. Вы можете сказать, что это не было сделано, потому что phi_j не дифференцируется в LHS.

2) Я бы предложил не использовать "А" в качестве дифференциального оператора. Вы также используете этот символ для матрицы глобальной системы, поэтому ваши обозначения противоречивы. Люди часто используют "D", так как это лучше соответствует идее, что оно используется для дифференциации.

Во-вторых, о вашей реализации:

3) Вы используете больше точек интеграции, чем необходимо. Ваши элементы используют функции линейной интерполяции, что означает, что вам нужна только одна точка интегрирования, расположенная в центре элемента, чтобы точно оценить интеграл. Изучите детали квадратуры Гаусса, чтобы понять почему. Также вы указали число точек интеграции, кратное количеству узлов. Вместо этого это должно быть кратно количеству элементов (в вашем случае n_vertex-1), потому что элементы - это домены, в которые вы интегрируетесь.

4) Вы создали свою систему, просто удалив два конечных узла из состава. Это не правильный способ указать граничные условия. Я бы предложил сначала построить полную систему и использовать один из типичных методов применения граничных условий Дирихле. Кроме того, подумайте, что будет означать ограничение двух узлов для дифференциального уравнения, которое вы пытаетесь решить. Какая функция существует, удовлетворяющая x' = 1, x(0) = 0, x(5) = 0? Вы перегрузили систему, пытаясь применить 2 граничных условия к дифференциальному уравнению первого порядка.

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

РЕДАКТИРОВАТЬ, чтобы учесть ваши изменения:

1) Если матрица A адресована с помощью A[row,col], то ваши индексы обратные. Вы должны интегрироваться с A [i, j] =...

2) Простой способ применить ограничение - заменить одну строку требуемым ограничением. Если вы хотите, чтобы x (0) = 0, например, установите A[0,j] = 0 для всех j, то установите A[0,0] = 1 и установите b[0] = 0. Это заменит один из уравнения с u_0 = 0. Сделайте это после интегрирования.