豆豉鲮鱼罐头 2022-02-22 04:06 采纳率: 50%
浏览 50
已结题

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

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

原要求是用托马斯算法写一个公式
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的哪一个图呢?)

  • 写回答

1条回答 默认 最新

  • joel_1993 2022-02-22 09:08
    关注

    你好,需要知道的是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])
    
    

    结果

    img

    斜率是 0.5244323436897312

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论
    1人已打赏

报告相同问题?

问题事件

  • 系统已结题 3月2日
  • 已采纳回答 2月22日
  • 赞助了问题酬金5元 2月22日
  • 创建了问题 2月22日

悬赏问题

  • ¥15 prism提示我reinstall prism 如何解决
  • ¥15 asp.core 权限控制怎么做,需要控制到每个方法
  • ¥20 while循环中OLED显示中断中的数据不正确
  • ¥15 这个视频里的stm32f4代码是怎么写的
  • ¥15 JNA调用DLL报堆栈溢出错误(0xC00000FD)
  • ¥15 请教SGeMs软件的使用
  • ¥15 自己用vb.net编写了一个dll文件,如何只给授权的用户使用这个dll文件进行打包编译,未授权用户不能进行打包编译操作?
  • ¥50 深度学习运行代码直接中断
  • ¥20 需要完整的共散射点成像代码
  • ¥15 编写vba代码实现数据录入工作