Forum - Python技术 - 线性方程组求解系数矩阵随结果更新的问题

线性方程组求解系数矩阵随结果更新的问题
tccyx 2015/05/05

要求解一个线性方程组AX=BY,右边是已知结果,要求解X。但因为边界条件的问题,系数矩阵B中最后一个值需要X中的结果来时时更新,但运算过程中B变成了nan,新手求解!谢谢!

01A_c=np.diagflat([-sigma for i in range(J-1)],-1)+\
02 np.diagflat([1.+sigma]+[1.+2.*sigma for i in range(J-2)]+[1.-sigma+lamb*sigma])+\
03 np.diagflat([-sigma for i in range(J-1)],1)
04
05B_c=np.diagflat([sigma for i in range(J-1)],-1)+\
06 np.diagflat([1.-sigma]+[1.-2.*sigma for i in range(J-2)]+[1.-sigma+lamb*sigma*C[-1]])+\
07 np.diagflat([sigma for i in range(J-1)],1)
08
09B=B_c.copy()
10
11C_record=[]
12C_record.append(C)
13
14for time in range(1,N):
15 C_new=np.linalg.solve(A_c,B.dot(C))
16 B[-1][-1]=float(1.-sigma+lamb*sigma*C_new[-1])
17 C=C_new
18 C_record.append(C)
RY 2015/05/12

请帖完整的可以运行的程序。

RY 2015/05/12

我测试了一下如下程序,C的值越来越大,最后超过了浮点数的范围:

01J=5
02sigma = 1.0
03lamb = 1.0
04C = np.arange(J)
05N = 20
06A_c=np.diagflat([-sigma for i in range(J-1)],-1)+\
07 np.diagflat([1.+sigma]+[1.+2.*sigma for i in range(J-2)]+[1.-sigma+lamb*sigma])+\
08 np.diagflat([-sigma for i in range(J-1)],1)
09
10B_c=np.diagflat([sigma for i in range(J-1)],-1)+\
11 np.diagflat([1.-sigma]+[1.-2.*sigma for i in range(J-2)]+[1.-sigma+lamb*sigma*C[-1]])+\
12 np.diagflat([sigma for i in range(J-1)],1)
13
14B=B_c.copy()
15
16C_record=[]
17C_record.append(C)
18
19for time in range(1,N):
20    C_new=np.linalg.solve(A_c,B.dot(C))
21    B[-1][-1]=float(1.-sigma+lamb*sigma*C_new[-1])
22    C=C_new
23    C_record.append(C)

下面是C_record的内容:

01[array([0, 1, 2, 3, 4]),
02 array([ 2.0952381 , 3.19047619, 6.47619048, 14.23809524, 33.23809524]),
03 array([ 58.79861786, 114.40675953, 279.04070835, 711.76298456,
04         1830.77205485]),
05 array([ 159882.6113512 , 319650.81594287, 798846.40391072,
06         2076341.26675356, 5428779.34657131]),
07 array([ 1.40341247e+12, 2.80682463e+12, 7.01706077e+12,
08          1.82443561e+13, 4.77160034e+13]),
09 array([ 1.08419856e+26, 2.16839712e+26, 5.42099280e+26,
10          1.40945813e+27, 3.68627511e+27]),
11 array([ 6.47077341e+53, 1.29415468e+54, 3.23538670e+54,
12          8.41200543e+54, 2.20006296e+55]),
13 array([ 2.30489382e+109, 4.60978764e+109, 1.15244691e+110,
14          2.99636197e+110, 7.83663899e+110]),
15 array([ 2.92442432e+220, 5.84884864e+220, 1.46221216e+221,
16          3.80175161e+221, 9.94304268e+221]),
17 array([ nan, nan, nan, inf, inf]),
18 array([ nan, nan, nan, nan, nan]),
19 array([ nan, nan, nan, nan, nan]),
20 array([ nan, nan, nan, nan, nan]),
21 array([ nan, nan, nan, nan, nan]),
22 array([ nan, nan, nan, nan, nan]),
23 array([ nan, nan, nan, nan, nan]),
24 array([ nan, nan, nan, nan, nan]),
25 array([ nan, nan, nan, nan, nan]),
26 array([ nan, nan, nan, nan, nan]),
27 array([ nan, nan, nan, nan, nan])]

如果把lamb设置为0.01就不会发散。你可以说说你要解决的是什么问题吗?

loading...