Метод Эйлера, схема python

Question

Редактировать: функция u [t] = [[0, 1], [-6,5]] * u с начальным условием u (0) = [[1], [1]].Точный ответ: u_t = [[1], [2]] * 2 * e ** (2 * t) - [[1], [3]] * e ** (3 * t)

Я работаю над проектом, и я заблудился относительно того, почему это не работает.Я планирую запустить схему метода Эйлера в следующем коде.

import numpy as np
from matplotlib import pyplot as plt

def Eulm():

    x0=0
    y0=1
    z0=1
    n=21
    xf=2
    y0=1
    z0=1
    w0=y0
    q0=z0
    w = [0] * (n+1)
    q = [0] * (n+1)
    w[0]=w0
    q[0]=q0
    x=np.linspace(x0,xf,n)
    y=2*np.e**(2*x)-np.e**(3*x)
    z=4*np.e**(2*x)-3*np.e**(3*x)
    L=[0]

    for i in range (1,n):
        deltax=(xf-x0)/(n-1)
        y0=y[0]
        z0=z[0]
        A=np.matrix([[1, -deltax], [6*deltax, 1-5*deltax]])
        G=np.linalg.inv(A)
        print(G)
        b=np.matrix([y[i-1], z[i-1]])
        b=b.transpose
        w[i]=G[0][0]*w[0][i]+G[0][1]*q[0][i-1]
        q[i]=G[1][0]*w[0][i-1]+G[1][1]*q[0][i-1]
        plt.plot(x,y)
        plt.xlabel('Time')
        plt.ylabel('Numerical Solutions')
        plt.title('Numerical Solutions with respect to Time')
        plt.show()

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

TypeError                                 Traceback (most recent call last)
<ipython-input-49-da43a543ec92> in <module>()
----> 1 Eulm()

    <ipython-input-47-022f3e9099b9> in Eulm()
         29         b=np.matrix([y[i-1], z[i-1]])
         30         b=b.transpose
    ---> 31         w[i]=G[0][0]*w[0][i]+G[0][1]*q[0][i-1]
         32         q[i]=G[1][0]*w[0][i-1]+G[1][1]*q[0][i-1]
         33         plt.plot(x,y)

TypeError: 'int' object is not subscriptable

Любая помощь будет принята с благодарностью.Спасибо!

Answer 1

Существует большое количество проблем с вашим кодом / теоретической реализацией метода Эйлера.Вот простая версия eulm, которая делает то, что должна:

import numpy as np
from matplotlib import pyplot as plt

def eulm(n=2001):
    x0=0
    y0=1
    z0=1
    xf=2
    wq = np.zeros((n, 2))
    wq[0] = y0,z0
    x = np.linspace(x0,xf,n)
    y = 2*np.exp(2*x) - np.exp(3*x)
    z = 4*np.exp(2*x) - 3*np.exp(3*x)
    A = np.array([[0, 1], [-6, 5]])
    dx = (x[1] - x[0])

    fig = plt.figure(figsize=(10,6))
    ax = fig.gca()
    ax.set_xlabel('Time')
    ax.set_ylabel('Numerical Solutions')
    ax.set_title('All Numerical Solutions with respect to Time,\n%d timesteps' % n)

    for i in range (1,n):
        wq[i] = wq[i - 1] + dx*(A @ wq[i - 1])
        if i == n-1:
            # add a legend in the final display
            ax.plot(x, wq[:, 0], c='C0', label='approx y')
            ax.plot(x, wq[:, 1], c='C1', label='approx z')
            ax.legend()
            fig.show()
        else:
            ax.plot(x, wq[:, 0], c='C0')
            ax.plot(x, wq[:, 1], c='C1')
            fig.show()

    fig = plt.figure(figsize=(10,6))
    ax = fig.gca()
    ax.set_xlabel('Time')
    ax.set_ylabel('Numerical Solutions')
    ax.set_title('Final Numerical Solutions with respect to Time,\n%d timesteps' % n)
    ax.plot(x, wq[:, 0], ':', c='k', lw=3, zorder=99, label='approx y')
    ax.plot(x, y, c='C2', lw=15, label='exact y')
    ax.plot(x, wq[:, 1], '--', c='k', lw=3, zorder=99, label='approx z')
    ax.plot(x, z, c='C3', lw=15, label='exact z')
    ax.legend()
    fig.show()

    return wq, x, y, z

Я свернул w и q в один массив из двух столбцов wq.Я избавился от обратной матрицы G (я не мог понять, что вы пытались с ней сделать, и в любом случае она не является частью стандартного метода Эйлера ).Я упростил вычисление значения wq на следующем шаге по времени, используя некую матричную нотацию от Numpy (оператор N 1012 * был недавно добавлен в Numpy и функционирует как оператор умножения матриц).

При достаточно большом n (и, следовательно, достаточно малом временном шаге dx), оценочные значения в wq[:, 0] и wq[:, 1] сходятся к точным значениям в y и z снебольшая ошибка (около ~1% при n=2001).Вот как выглядят графики, которые eulm() генерирует при запуске: