前两天我问了个关于求解方程的问题，RY给了我一个很好的重复设置初值的方法，很有启发作用，改进后我试着编到自己的方程中，可发现导入的数值是元组类型，我需要的是可变的数组，自己试着修改，没改成功。把编写的程序发上来，希望大家能帮忙改正。还有个问题就是如何设置迭代方程的个数，比如我要求解20个或更多的方程，方程没有变化，只是其中几个数值有变化，可不可以它设置成和N有关循环的方程。

# -*- coding: cp936 -*-
import numpy as np
from scipy.optimize import leastsq
from scipy.optimize import fsolve
import pylab as pl
from math import exp,log
import random

def func(H):
    x=[i for i in xrange(N1)]    #给各个变量定义一个长度为N1的一个一维数组
    ESx=np.array(x,dtype=np.float32)

    for i in range(N1):
        ESx[i]= H[i]# 以下是迭代方程组，求非线性方程
        return [ V - (ESx[0]*L[0]+ESx[1]*L[1]+ESx[2]*L[2]+ESx[3]*L[3]+ESx[4]*L[4] ),
            (S[0] + k*(ESx[0]*dt)**n + m*log(ESx[0]))*AS[0]/exp(ESx[0]*dt)-(S[1]\
             + k*(ESx[1]*dt)**n + m*log(ESx[1]))*AS[1]/exp(ESx[1]*dt),
            (S[1] + k*(ESx[1]*dt)**n + m*log(ESx[1]))*AS[1]/exp(ESx[1]*dt)-(S[2]\
             + k*(ESx[2]*dt)**n + m*log(ESx[2]))*AS[2]/exp(ESx[2]*dt),     
            (S[2] + k*(ESx[2]*dt)**n + m*log(ESx[2]))*AS[2]/exp(ESx[2]*dt)-(S[3]\
             + k*(ESx[3]*dt)**n + m*log(ESx[3]))*AS[3]/exp(ESx[3]*dt),
            (S[3] + k*(ESx[3]*dt)**n + m*log(ESx[3]))*AS[3]/exp(ESx[3]*dt)-(S[4]\
             + k*(ESx[4]*dt)**n + m*log(ESx[4]))*AS[4]/exp(ESx[4]*dt) ]
                   

A0=10  #原始面积
L0=2  #原始长度
L_total=10  # 总长度
V=0.1  # 拉伸速度
k=0.5  # 塑性参数
m=0.3  # 应变敏感指数
n=0.2  # 应变硬化指数
N1=5   # 有限单元数
M1=2   # 颈缩单元数
dt=0.25  # 时间增量
fa0=0.005  # 颈缩处面积比
fra_cond=2 # 断裂条件
S=np.random.uniform(10,20,N1) # 创建一个一维数组，存放N1个随机的初始应力增量
c=[]

x=[i for i in xrange(N1)]    #给各个变量定义一个长度为N1的一个一维数组
f=np.array(x,dtype=np.float32)
L=np.array(x,dtype=np.float32)
E=np.array(x,dtype=np.float32)
ES=np.array(x,dtype=np.float32)
AS=np.array(x,dtype=np.float32)
A=np.array(x,dtype=np.float32)

b=np.arange(N1)    #所有单元拍成列表放在b中

random.shuffle(b)  #洗牌 防止不同次随机出现的值相同  
for i in range(M1):
    c.append(b[i])      #随机颈缩出现的位置
    a=c[i]
    f[a]=(1-fa0)-fa0*c[i] #f即为保存随机颈缩和零值的列表
##print f
#以上是在进行随机设置颈缩并复制随机的初值

for i in range(N1):  # 原始值
    L[i]=L0
    E[i]=0          # 给每个单元赋予初始应变量
    ES[i]=V/(L0*N1)  #给每个单元赋予初始应变率
    AS[i]=A0
    for j in range(M1):
        if c[j]==i:
            AS[i]=f[j]*A0 # 赋予随即颈缩

ES_total=sum(ES)
while (max(ES)/ES_total<fra_cond):   # 断裂条件
    for j in range(N1):
        # 以下都是进行跟时间有关系的迭代运算，求解不同时刻的值
        E[j]=ES[j]*dt+E[j] # 真实应变
        L[j]=L0*np.exp(E[j])   # 真实长度
        A[j]=AS[j]/np.exp(E[j])   # 真实面积
       
   
    H = ES
    ES=fsolve(func,H)    # fn返回真实的应变率
    "fsolve是求解非线性方程的内嵌函数，返回值类型和H一样"
    print ES
   
    L_total=sum(L)             #  total tensile length
    E_total=log(sum(L)/(L0*N1))  # total tensile strain
    ES_total=V/(sum(L))         # total strain rate
    p=(S[0] + k*(ESx[0]*dt)**n + m*log(ESx[0]))*AS[0]/exp(ESx[0]*dt)     # load
    print p
##    d_total=p/((l0*N1*a0)/l_total)  # total stress
    time=time+dt    # tensile time
print time
  

已在 http://hyry.dip.jp/tech/forum/thread.html/54 回复，

你需要用NumPy的矢量运算功能将程序简化。

谢谢你的宝贵意见，我对Numpy的的函数知道的太少了， 还需要学习啊。

现在这个问题还是出在fsolve函数上，我用你提供的建议，用tmp[:-1] - tmp[1:] 来代替我的那些方程，可效果并不好，提示类型错误“there is a mismatch between the input and output shape of the 'func' argument 'func'”，“我的个人理解是把tmp[:-1] - tmp[1:] 当成了一个方程，而不是N-1个方程来求解。”如果用我的把所有的方程列出来求解，就没有这个问题，但是迭代效果不好，结果离真值太远，这两个问题都应该怎么解决啊。卡在这个函数上N多天了，麻烦你了。

import numpy as np
from scipy.optimize import fsolve
def func(H):
    tmp =(H*0.5)**0.2 + 1000*np.log(H)*AS/np.exp(H*0.5)
    return [tmp[:-1]-tmp[1:],
            10-sum(H)]
##    return[(H[0]*0.5)**0.2 + 1000*np.log(H[0])*AS[0]/np.exp(H[0]*0.5)-(H[1]*0.5)**0.2 + 1000*np.log(H[1])*AS[1]/np.exp(H[1]*0.5),
##           (H[1]*0.5)**0.2 + 1000*np.log(H[1])*AS[1]/np.exp(H[1]*0.5)-(H[2]*0.5)**0.2 + 1000*np.log(H[2])*AS[2]/np.exp(H[2]*0.5),
##            10-sum(H)]
   
ES=np.array([1,2,3],dtype=float)
AS=np.array([10,11,12],dtype=float)
ES=fsolve(func,ES)
print ES

运算完毕之后还是要把各个方程的结果还原成一个列表：

return list(tmp[:-1]-tmp[1:]) + [10-sum(H)]

至于迭代效果不好则只能多用一些初始值，然后选取误差最小的结果了。

你给我得return list(tmp[:-1]-tmp[1:]) + [10-sum(H)]不但不出错误了，还迭代得到好的结果了。但应用了更复杂的方程，就比如我得那个方程，效果就不好了。现在问题已经解决一半，目前只剩下初值的设置了，呵呵。

你之前给过我一个关于反复设置初值设置的建议，用下面的这个函数得到了一个2维的坐标作为初值不断地尝试，是个很好的想法，不过不太适合我的这个输入数值是N维的运算，而且我想要的返回值是一个N维的数组，不是个元组。

我想尽可能的找到一个有规律的方法来设置初值，就比如你给出的这种尝试着作为初值来解的方法。不知道你还有没有好的建议

tmp = np.arange(0.1, 1.0, 0.1)

for p0 in itertools.product(tmp, tmp)

N元方程组的初值问题的确比较麻烦，因为循环次数会很多，可能需要寻找根据方程的公式推测出大概的初值范围的方法。

实在不行也可以试试随机数，np.random.normal(1.0, 2.0, 10)，这个产生10个均值为1，标准差为2.0的10个随机数。

如果你要在N维网格上设置初始值的话：

tmp = np.arange(0.1, 1.0, 0.2)

for p0 in itertools.product(*[tmp]*n):

    p0 = np.array(p0)