1

我在以下代码中遇到错误，关于矩阵，有人可以解释发生了什么吗？我是 python 新手，无法真正理解错误。该代码是使用凯勒盒求解边界层，也可以称为矩阵中的矩阵。这是整个代码

import numpy as npy

blt = int(raw_input("Input the boundary layer thickness = "))
deleta = float(raw_input("Input the step size of boundary layer thickness = "))
np = int((blt/deleta) + 1)
stop = 1
k=1
l=2
g=2
eselon = 0.00001

def eta(j,k):
    if j == 1 and k == 1:
        return 0
    else:
        return eta(j-1,k) + deleta;
deta = deleta
def etanp():
    return eta(np,k)       
def f(j,k):
    return -eta(np,1)*etab2 + etanp()*etab1 + eta(j,1)
def u(j,k):
    return  -3*etab1 + 2*etab +1
def v(j,k):
    return (-6/etanp()*etab + (2/etanp()))
def fb(j,k):
    return 0.5 * (f(j,k) + f(j-1,k))
def ub(j,k):
    return 0.5 * (u(j,k) + u(j-1,k))
def vb(j,k):
    return 0.5*(v(j,k) + v(j-1,k))
def fvb(j,k):
    return fb(j,k)*vb(j,k)
def uub(j,k):
    return ub(j,k) * ub(j,k)

def a1(j,k):
    return 1 + 0.5*deta *fb(j,k)
def a2(j,k):
    return -1 + 0.5*deta *fb(j,k)
def a3(j,k):
    return 0.5 * deta * vb(j,k)
def a4(j,k):
    return a3(j,k)
def a5(j,k):
    return -1 * deta * ub(j,k)
def a6(j,k):
    return a5(j,k)

def r1(j,k):
    return f(j-1,k) - f(j,k) + deta * ub(j,k)
def r2(j,k):
    return u(j-1,k) - u(j,k) + deta * vb(j,k)
def r3(j,k):
    return v(j-1,k)-v(j,k) - deta *((fvb(j,k)*uub(j,k)))-deta

def AJ(j,k):
    if j == 2:
        return npy.matrix([[0,1,0],[-0.5*deta,0,-0.5*deta],[a2(2,k), a3(2,k), a1(2,k)]])
    else:
        return npy.matrix([[-0.5*deta,0,0],[-1,0,-0.5*deta],[a6(j,k),a3(j,k),a1(j,k)]])
def BJ(j,k):
    return npy.matrix([[0,-1,0],[0,0,-0.5*deta],[0,a4(j,k),a2(j,k)]])
def CJ(j,k):
    return npy.matrix([[-0.5*deta,0,0],[1,0,0],[a5(j,k),0,0]])

def alfa(j,k):
    return AJ(j,k) - (BJ(j,k)*gamma(j,k))
def gamma(j,k):
    return npy.matrix.I((alfa(g,k))*CJ(g,k))
def rr(j,k):
    return npy.matrix([[r1(j,k)],[r2(j,k)],[r3(j,k)]])   
def ww(j,k):
    if j == 2:
        return npy.matrix.I(AJ(2,k)*rr(2,k))
    else:
        return npy.matrix.I((alfa(j,k))*(rr(j,k)-(BJ(j,k)*ww(j-1,k))))

def dell(j,k):
     if j == np:
        return ww(np,k)
     else:   
        return ww(j,k) - (gamma(j,k)*dell(j+1,k))
def delf(j,k):
    if j == 1:
        return 0
    elif j == 2:
        return dell(2,k)[2,1]
    else:
        return dell(j,k)
def delu(j,k):
    if j == 1 or j == np:
        return 0
    elif j == np-1:
        return dell(j,k)[1,1]
def delv(j,k):
    if j == 1:
        return dell(2,k)[1,1]
    elif j == 2:
        return dell(2,k)[3,1]
    else:
        return dell(j,k)[3,1]

def ffinal(j,l):
    return f(j,k) + delf(j,k)
def ufinal(j,l):
    return u(j,k) + delu(j,k)
def vfinal(j,l):
    return v(j,k) + delv(j,k)

# Beginning of calculation for Keller-Box

while stop > eselon:
    eta(1,1)
    for j in range (2,np):
        eta(j,k)  

# Initial condition
    etab = eta(j,k) / eta(np,k)
    etab1 = etab**2
    etab2 = etab**3 
    for j in range (1,np):
        deta
        f(j,1)
        u(j,1)
        v(j,1)

# Current value of Central Differentiation
    for j in range (2,np):
        fb(j,k)
        ub(j,k)
        vb(j,k)
        fvb(j,k)
        uub(j,k)
        a1(j,k)
        a2(j,k)
        a3(j,k)
        r1(j,k)
        r2(j,k)
        r3(j,k)
# Matrices Value for A1, Aj, Bj, and CJ
        CJ(j,k)
        AJ(j,k)
        BJ(j,k)
# Recursion: Forward Sweeping
    for j in range (3,np):
        alfa(j,k)
        gamma(j,k)
    for j in range(2,np):
        rr(j,k)
    for j in range(3,np):
        ww(j,k)

# Recursion: Backward Sweeping
    for j in range (np-1,2,-1):
        dell(j,k)

    for j in range (np,3,-1):
        delu(j-1,k)
        delf(j,k)
        delv(j,k)

# Newton's Method
    for j in range (1,np):
        ffinal(j,l)
        ufinal(j,l)
        vfinal(j,l)

# Check the convergence of iteration
    stop = npy.abs(delv(1,k))
    kmax = k
    k =+ 1

cfrex = vfinal(1,kmax)

print cfrex
4

2 回答 2

1

你应该用括号调用 numpy.matrix ：

>>> numpy.matrix([[1,2],[1,2]])
matrix([[1, 2],
        [1, 2]])

试试numpy 教程。

于 2013-04-12T11:28:01.900 回答
1

更新：您的函数alfa和gamma，定义如下：

def alfa(j,k):
    print 'alfa({},{}) called'.format(j,k)
    return AJ(j,k) - (BJ(j,k)*gamma(j,k))

def gamma(j,k):
    print 'gamma({},{}) called'.format(j,k)
    return npy.matrix.I((alfa(g,k))*CJ(g,k))

无论它们作为参数给出什么，都将永远运行，因为每个都只会调用另一个并且它们无法终止。也就是说，无论输入如何，每个函数都会调用另一个函数。

我可以看到您的代码中有一些错误。

1）这是@graphite提到的一个，在npy.matrix其参数周围没有括号的情况下被调用：

def A1(j,k):
    return npy.matrix[[0,1,0],[-0.5*deta(2,k),0,-0.5*deta(2,k)],[a2(2,k), a3(2,k), a1(2,k)]] 
# needs parentheses:
    return npy.matrix([[0,1,0],[-0.5*deta(2,k),0,-0.5*deta(2,k)],[a2(2,k), a3(2,k), a1(2,k)]])

在您使用的任何地方修复此问题npy.matrix

2）在您的定义中r，您缺少我认为应该是的内容k：

def r1(j,k):
    return f(j-1,k) - f(j,k) + deta(j,) * ub(j,k)
#should be:
    return f(j-1,k) - f(j,k) + deta(j,k) * ub(j,k)

3) 这里np1没有定义。也许应该是np？

# Recursion: Backward Sweeping
    for j in range (np1,2,-1):
        dell(j,k)

但总的来说，如果你定义的函数更少，尤其是那些没有参数的函数，那么这段代码可以变得更具可读性（因此更容易找到错误），例如：

def etab():
    return eta(j,k) / eta(np,k)
def etab1():
    return etab()*etab()
def etab2():
    return etab()*etab()*etab()

可以替换为：

etab = eta(j, k) / eta(np, k)
etab1 = etab**2
etab2 = etab**3

另一个例子：

def deta(j,k):
    return deleta

甚至不使用jor k，因此您可以将整个内容更改为：

deta = deleta

在下面的任何地方，而不是写作，deta(j, k)或者deta(2, k)你可以只写deta. 当然，如果您计划在未来某个时候更改您的代码，以便deta实际更改为不同的jor k，那么您可以保持原样。

我希望这有帮助！

于 2013-04-12T22:14:12.083 回答