Matlab转成python码 一直在报错 请问能帮我改一下吗？

Question

问题遇到的现象和发生背景

原要求是用托马斯算法写一个公式
4u - d2u/dx2 = 0,
boundary是 u(0) =0, u(1) = 10
要求算N = 10 20 40 80分别的analytical solution (ua) 和 central difference solution (ui)
然后求两个解的error
err=sqrt(sum(diag(ui-ua).^2));

问题相关代码，请勿粘贴截图

这是我原先的matlab码

N = 10; %10, 20, 40,  80
Num = N+2;
h = 1.0/(Num-1);
x0=ones(Num,1);
x1=ones(Num-1,1);
A=4*eye(Num)-(diag(-2*x0)+diag(x1,1)+diag(x1,-1))/(h*h);
A(1,:)=0;
A(Num,:)=0;
A(1,1)=1;
A(Num,Num)=1;
b=zeros(Num,1);
b(1)=0;
b(Num)=10;
x=linspace(0,1,Num);
ui=inv(A)*b;
AA=(-10/((exp(-2)-exp(2))))*exp(2*x)+(10/((exp(-2)-exp(2))))*exp(-2*x);
ua=AA
plot(x,ui,x,ua,'Marker','*')
single_err=diag(ui-ua);
sqr_err=single_err.^2;
err=sqrt(sum(sqr_err));
err;
%now I'll try to convert everything to horizontal format and export ua, ui,
%and x onto the same sheet. %
sheet = [x;ui';ua];

最后可以对每一个N（10，20，40，80）画出一个ua ui的线图，同时给出对应的error值。我想把这个码改写成python形式。

运行结果及报错内容

ValueError Traceback (most recent call last)
in
----> 1 ThomasAGiven(10)

in ThomasAGiven(N)
7 x0 = np.ones((Num,1))
8 x1 = np.ones((Num-1,1))
----> 9 A = 4np.eye(Num)-(np.diag(-2x0)+np.diag(x1,1)+np.diag(x1,-1))/(h*h)
10 A[1,:] = 0
11 A[Num,:] = 0

ValueError: operands could not be broadcast together with shapes (12,12) (0,)

我的解答思路和尝试过的方法

以下是我当前改的，但一直说我operand有问题，向量/矩阵格式不对

def ThomasAGiven(N):
    import numpy as np
    import matplotlib.pyplot as plt 
    Num = N + 2
    h = 1/(Num-1)

    x0 = np.ones((Num,1))
    x1 = np.ones((Num-1,1))
    A = 4*np.eye(Num)-(np.diag(-2*x0)+np.diag(x1,1)+np.diag(x1,-1))/(h*h)
    A[1,:] = 0
    A[Num,:] = 0
    A[1,1] = 1
    A[Num,Num] = 1
    b = np.zeros[Num,1]
    b[1] = 0
    b[Num] = 10
    x = np.linspace(0,1,Num)
    ui = np.linalg.inv(A)*b

    AA = (-10/((np.exp(-2)-np.exp(2))))*np.exp(2*x)+(10/((np.exp(-2)-np.exp(2))))*np.exp(-2*x)
    ua = AA
    single_err = np.diag(ui-ua)
    sqr_err = single_err**2
    err = np.sqrt(np.sum(sqr_err))
    return err

我想要达到的结果

改成python码，同时把每一个N对应的error值写入注明斜率的loglog图。
（log每一个N值，同时log每一个对应的error值，然后插入图表，这个线状图应该用matplotlib的哪一个图呢？）

Answer 1

你好，需要知道的是python的下标从0开始的

def ThomasAGiven(N):
    import numpy as np
    Num = N + 2
    h = 1/(Num-1)
 
    x0 = np.ones(Num)
    x1 = np.ones(Num-1)
    A = 4*np.eye(Num)-(np.diag(-2*x0)+np.diag(x1,1)+np.diag(x1,-1))/(h*h)
    A[0,:] = 0
    A[Num-1,:] = 0
    A[0,0] = 1
    A[Num-1,Num-1] = 1
    b = np.zeros((Num,1))
    b[0] = 0
    b[Num-1] = 10
    x = np.linspace(0,1,Num)
    ui = np.linalg.inv(A)*b
    AA = (-10/((np.exp(-2)-np.exp(2))))*np.exp(2*x)+(10/((np.exp(-2)-np.exp(2))))*np.exp(-2*x)
    ua = AA
    single_err = np.diag(ui-ua)
    sqr_err = single_err**2
    err = np.sqrt(np.sum(sqr_err))
    return err
import numpy as np
import matplotlib.pyplot as plt 
N = np.array([10,20,40,80])
err = np.zeros(4)
for i in range(len(N)):
    err[i] = ThomasAGiven(N[i])
plt.plot(np.log(N),np.log(err))
plt.show()
p = np.polyfit(np.log(N),np.log(err),1)%线性拟合
print('斜率是', p[0])

结果

斜率是 0.5244323436897312