这是非语法错误代码，但我似乎无法修复递归错误。在这里需要一些帮助。基于 matlab 的算法，我已经阅读了 matlab 上的教程，但我似乎可以弄清楚我错过了哪一部分。
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):
if j == 2:
return AJ(2,k)
else:
return AJ(j,k) - (BJ(j,k)*gamma(j-1,k))
def gamma(j,k):
if j == 2:
return npy.matrix.I((alfa(2,k))*CJ(2,k))
else:
return npy.matrix.I((alfa(j,k))*CJ(j,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)[1,0]
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)[0,0]
def delv(j,k):
if j == 1:
return dell(2,k)[0,0]
elif j == 2:
return dell(2,k)[2,0]
else:
return dell(j,k)[2,0]
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
这是我从 mathlab 使用的参考资料
*******************************************************************
Input
*******************************************************************
blt = input ('Input the boundary layer thickness = ');
deleta=0.1; %input('Input the step size of boundary layer thickness=');
np = (blt / deleta)+ 1;
pr = 7; %input ('Input the prandtl number = ');
K = 0; %input ('Input the material parameter K = ');
lambda = 1; %input ('Input the mixed convection parameter = ');
stop = 1.0; k = 1; eselon = 0.00001;
while stop > eselon
eta(1,1) = 0.0;
for j = 2:np
eta(j,1) = eta(j-1,1) + deleta;
end
*******************************************************************
Initial Condition for f, u, v, h, p, s, t
*******************************************************************
etanpq = eta(np,1)/3;
etau15 = 1/eta(np,1);
etau16 = 2/eta(np,1);
etanp = eta(np,1);
for j =1:np
deta(j,k)= deleta;
etab = eta(j,1)/eta(np,1);
etab1 = etab^2;
etab2 = etab^3;
etau152 = etau15^2;
etau162 = etau16^2;
f(j,1) = -etanpq * etab2 + etanp * etab1;
u(j,1) = -etab1 + 2 * etab;
v(j,1) = -etau16 * etab + etau16;
h(j,1) = etau15 * etab - etau15;
p(j,1) = etau152;
s(j,1) = -eta(j,1) + eta(j,1) * etab;
t(j,1) = -1 + 2 * etab;
end
*******************************************************************
Current Central Differention Value
*******************************************************************
for j = 2:np
fb(j,k) = 0.5*(f(j,k) + f(j-1,k));
ub(j,k) = 0.5*(u(j,k) + u(j-1,k));
vb(j,k) = 0.5*(v(j,k) + v(j-1,k));
hb(j,k) = 0.5*(h(j,k) + h(j-1,k));
pb(j,k) = 0.5*(p(j,k) + p(j-1,k));
sb(j,k) = 0.5*(s(j,k) + s(j-1,k));
tb(j,k) = 0.5*(t(j,k) + t(j-1,k));
fvb(j,k) = fb(j,k) * vb(j,k);
uub(j,k) = ub(j,k) ^ 2;
pfb(j,k) = pb(j,k) * fb(j,k);
hub(j,k) = hb(j,k) * ub(j,k);
tfb(j,k) = tb(j,k) * fb(j,k);
sub(j,k) = sb(j,k) * ub(j,k);
*******************************************************************
Momentum Differential Equation
*******************************************************************
a1(j,k) = (1.0 + K) + 0.5 * deta(j,k) * fb(j,k);
a2(j,k) = -(1.0 + K) + 0.5 * deta(j,k) * fb(j,k);
a3(j,k) = 0.5 * deta(j,k) * vb(j,k);
a4(j,k) = a3(j,k);
a5(j,k) = -1 * deta(j,k) * ub(j,k);
a6(j,k) = a5(j,k);
a7(j,k) = 0.5 * K * deta(j,k);
a8(j,k) = a7(j,k);
a9(j,k) = 0.5 * lambda * deta(j,k);
a10(j,k) = a9(j,k);
*******************************************************************
Angel Differential
*******************************************************************
b1(j,k) = (1 + K/2) + 0.5 * deta(j,k) * fb(j,k);
b2(j,k) = -(1 + K/2) + 0.5 * deta(j,k) * fb(j,k);
b3(j,k) = 0.5 * deta(j,k) * pb(j,k);
b4(j,k) = b3(j,k);
b5(j,k) = -0.5 * deta(j,k) * hb(j,k);
b6(j,k) = b5(j,k);
b7(j,k) = -0.5 * deta(j,k) * ub(j,k) - K * deta(j,k);
b8(j,k) = b7(j,k);
b9(j,k) = -0.5 * K * deta(j,k);
b10(j,k) = b9(j,k);
*******************************************************************
Energy Differential
*******************************************************************
c1(j,k) = 1/pr + 0.5 * deta(j,k) * fb(j,k);
c2(j,k) = -1/pr + 0.5 * deta(j,k) * fb(j,k);
c3(j,k) = 0.5 * deta(j,k) * tb(j,k);
c4(j,k) = c3(j,k);
c5(j,k) = -0.5 * deta(j,k) * sb(j,k);
c6(j,k) = c5(j,k);
c7(j,k) = -0.5 * deta(j,k) * ub(j,k);
c8(j,k) = c7(j,k);
*******************************************************************
Definition value of rj-1/2
*******************************************************************
r1(j,k) = f(j-1,k) - f(j,k) + deta(j,k) * ub(j,k);
r2(j,k) = u(j-1,k) - u(j,k) + deta(j,k) * vb(j,k);
r3(j,k) = h(j-1,k) - h(j,k) + deta(j,k) * pb(j,k);
r4(j,k) = s(j-1,k) - s(j,k) + deta(j,k) * tb(j,k);
r5(j,k) = (1.0 + K) * (v(j-1,k) - v(j,k)) - deta(j,k) * fvb(j,k) -...
deta(j,k)+ deta(j,k) * uub(j,k) - K * deta(j,k)...
* pb(j,k) - lambda * deta(j,k) * sb(j,k);
r6(j,k) = (1 + K/2) * (p(j-1,k) - p(j,k)) - deta(j,k) * pfb(j,k) + ...
deta(j,k) * hub(j,k) + 2 * K * deta(j,k) * hb(j,k) + ...
K * deta(j,k) * vb(j,k);
r7(j,k) = 1/pr * (t(j-1,k) - t(j,k)) - deta(j,k) * tfb(j,k) +...
deta(j,k) * sub(j,k);
end
*******************************************************************
Matrices Value A1, Aj, Bj, Cj
*******************************************************************
a{2,k} = [0 0 0 1 0 0 0 ;...
-0.5 * deta(2,k) 0 0 0 -0.5 * deta(2,k) 0 0;...
0 -0.5 * deta(2,k) 0 0 0 -0.5 * deta(2,k) 0;...
0 0 -1 0 0 0 -0.5 * deta(2,k);...
a2(2,k) a8(2,k) a10(2,k) a3(2,k) a1(2,k) a7(2,k) 0;...
b10(2,k) b2(2,k) 0 b3(2,k) b9(2,k) b1(2,k) 0;...
0 0 c8(2,k) c3(2,k) 0 0 c1(2,k)];
for j = 3:np
a{j,k} = [-0.5 * deta(j,k) 0 0 1 0 0 0 ;...
-1 0 0 0 -0.5 * deta(j,k) 0 0 ;...
0 -1 0 0 0 -0.5 * deta(j,k) 0 ;...
0 0 -1 0 0 0 -0.5 * deta(j,k);...
a6(j,k) 0 a10(j,k) a3(j,k) a1(j,k) a7(j,k) 0;...
b6(j,k) b8(j,k) 0 b3(j,k) b9(j,k) b1(j,k) 0;...
c6(j,k) 0 c8(j,k) c3(j,k) 0 0 c1(j,k)];
b{j,k} = [0 0 0 -1 0 0 0 ;...
0 0 0 0 -0.5 * deta(j,k) 0 0;...
0 0 0 0 0 -0.5 * deta(j,k) 0;...
0 0 0 0 0 0 -0.5 * deta(j,k);...
0 0 0 a4(j,k) a2(j,k) a8(j,k) 0;...
0 0 0 b4(j,k) b10(j,k) b2(j,k) 0;...
0 0 0 c4(j,k) 0 0 c2(j,k)];
end
for j = 2:np
c{j,k} = [-0.5 * deta(j,k) 0 0 0 0 0 0 ;...
1 0 0 0 0 0 0;...
0 1 0 0 0 0 0;...
0 0 1 0 0 0 0;...
a5(j,k) 0 a9(j,k) 0 0 0 0;...
b5(j,k) b7(j,k) 0 0 0 0 0;...
c5(j,k) 0 c7(j,k) 0 0 0 0];
end
*******************************************************************
Recursion of block Elimination
*******************************************************************
Forward Sweeping
*******************************************************************
alfa{2,k} = a{2,k};
gamma{2,k} = inv(alfa{2,k}) * c{2,k};
for j = 3:np
alfa{j,k} = a{j,k} - (b{j,k} * gamma{j-1,k});
gamma{j,k} = inv(alfa{j,k}) * c{j,k};
end
for j = 2:np
rr{j,k} = [r1(j,k); r2(j,k); r3(j,k); r4(j,k); r5(j,k); r6(j,k);...
r7(j,k)];
end
ww{2,k} = inv(alfa{2,k}) * rr{2,k};
for j = 3:np
ww{j,k} = inv(alfa{j,k}) * (rr{j,k} - (b{j,k} * ww{j-1,k}));
end
*******************************************************************
Backward Sweeping
*******************************************************************
delf(1,k) = 0.0;
delu(1,k) = 0.0;
delh(1,k) = 0.0;
delt(1,k) = 0.0;
delu(np,k) = 0.0;
delh(np,k) = 0.0;
dels(np,k) = 0.0;
dell{np,k} = ww{np,k};
for j = np-1:-1:2
dell{j,k} = ww{j,k} -(gamma{j,k} * dell{j+1,k});
end
delv(1,k) = dell{2,k}(1,1);
delp(1,k) = dell{2,k}(2,1);
dels(1,k) = dell{2,k}(3,1);
delf(2,k) = dell{2,k}(4,1);
delv(2,k) = dell{2,k}(5,1);
delp(2,k) = dell{2,k}(6,1);
delt(2,k) = dell{2,k}(7,1);
for j = np:-1:3
delu(j-1,k) = dell{j,k}(1,1);
delh(j-1,k) = dell{j,k}(2,1);
dels(j-1,k) = dell{j,k}(3,1);
delf(j,k) = dell{j,k}(4,1);
delv(j,k) = dell{j,k}(5,1);
delp(j,k) = dell{j,k}(6,1);
delt(j,k) = dell{j,k}(7,1);
end
*******************************************************************
Newton method
*******************************************************************
for j = 1:np
f(j,k+1) = f(j,k) + delf(j,k);
u(j,k+1) = u(j,k) + delu(j,k);
v(j,k+1) = v(j,k) + delv(j,k);
h(j,k+1) = h(j,k) + delh(j,k);
p(j,k+1) = p(j,k) + delp(j,k);
s(j,k+1) = s(j,k) + dels(j,k);
t(j,k+1) = t(j,k) + delt(j,k);
h(j,k+1) = -0.5 * v(j,k+1);
end
*******************************************************************
Convergence Check
*******************************************************************
stop = abs(delv(1,k));
kmax = k;
k = k + 1;
end
*******************************************************************
Skin Friction and Nusselt Number
*******************************************************************
cfrex = v(1,kmax)
nuxrex = 1/s(1,kmax)
nuxrex2 = s(1,kmax)
但是，我的研究案例仅涵盖 f、u 和 v，因此，初始条件有一些变化。和许多其他部分。