欧拉方法方案python

Question

编辑：函数为u[t]=[[0，1]，[-6,5]]*u，初始条件为u（0）=[[1]，[1]]。[t*2]，[t*2]，[t*2]，[t*2]

我正在做一个项目，我不知道为什么这不起作用。我计划在下面的代码上运行一个Euler方法方案。在

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()

我得到了下面的错误，我不确定为什么我会得到这个错误，考虑到它应该在列表中迭代。在

^{pr2}$

任何帮助都将不胜感激。谢谢您！在

Answer 1

Euler方法的代码/理论实现存在大量问题。下面是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（我不知道你想用它做什么，无论如何它不是{a1}的一部分）。我使用Numpy的一些矩阵表示法简化了下一个时间步中wq值的计算（Numpy最近添加了@运算符，并用作矩阵乘法运算符）。在

对于足够大的n（因此时间步长足够小dx），在wq[:, 0]和{}中的估计值收敛到y和{}中的精确值，误差很小（大约在n=2001时~1%）。运行eulm()时生成的绘图如下：