Метод Криво-Никольсона - PullRequest
Новая
       21

Метод Криво-Никольсона

0 голосов
/ 09 марта 2020

Мне нужно написать следующий псевдокод в Python код: введите описание изображения здесь

И вот мой код: 

import math

def f(x):
    v = math.sin(math.pi/2*x)
    return v

def zero_matrix(i, j, start_from=1):
    matrix = [[-1] * (j+start_from)]
    for x in range(i):
        line = [-1] * start_from + [0] * j
        matrix.append(line)
    return matrix

def zero_vector(i, start_from=1):
    return [-1] * start_from + [0] * i

def algo_12_3(l, T, alpha, m, N):
    """Crank-Nicolson Method.

    Args:
        l: endpoint
        T: maximum time
        alpha: constant
        m:
        N:
    Return:
        approximations wi,j to u(xi, tj) for each i=i,...,m-1 and j=1,...,N.
    """
    w = zero_matrix(m, N)
    z = zero_vector(m-1)
    ll = zero_vector(m-1)
    u = zero_vector(m-1)

    h = l / m
    k = T / N
    lambda_ = alpha * alpha * k / (h * h)
    for i in range(1, m):
        w[i][0] = f(i * h)
    ll[1] = 1 + lambda_
    u[1] = -lambda_/(2*ll[1])
    for i in range(2, m-1):
        ll[i] = 1 + lambda_ + lambda_ * u[i-1]/2
        u[i] = -lambda_/(2 * ll[i])
        ll[m-1] = 1 + lambda_ + lambda_ * u[m-2]/2
    for j in range(1, N + 1):
        t = j * k
        z[1] = ((1 - lambda_) * w[1][j-1] + lambda_ / 2 * w[2][j-1]) / ll[1]
        for i in range(2, m):
            z[i] = ((1 - lambda_) * w[i][j-1] + lambda_ / 2 * (w[i+1][j-1] + w[i-1][j-1] + z[i-1])) / ll[i]
        w[m-1][j] = z[m-1]
        for i in range(m-2, 0, -1):
            w[i][j] = z[i] - u[i] * w[i+1][j]
        print('t={}: '.format(t))
        for i in range(1, m):
            print('({}, {})'.format(i*h, w[i][j]))

algo_12_3(2, 0.1, 1, 4, 2)

Мои выводы : 

t=0.05: 
(0.5, 0.6282614874855517)
(1.0, 0.8822836003342388)
(1.5, 0.5449281541522184)
t=0.1: 
(0.5, 0.55687611895338)
(1.0, 0.7741379272219071)
(1.5, 0.5013205633978245)

Однако правильные выходные данные должны быть: 

t=0.05: 
(0.5, 0.62884835)
(1.0, 0.88932586)
(1.5, 0.62884835)
t=0.1: 
(0.5, 0.59404972)
(1.0, 0.84011317)
(1.5, 0.59404972)

Я не знаю, с диапазоном или с формированием исходной матрицы. Может ли кто-нибудь помочь мне с этим? Большое спасибо!! Цени это!

...