求解线性方程组

麦田守望者

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	SUBROUTINE SLBSI(A,B,D,N,MS,NX,MX)

02	DIMENSION A(NX,MX),B(NX),D(MX)

03

N1=N-1

04

DO K=1,N1

05

C=A(K,1)

06

K1=K+1

07	IF((ABS(C)-0.000001).LE.0)THEN

08	WRITE(6,2)K

09	2  FORMAT('*****SINGULARITY IN ROW',I5)

10

STOP

11

ELSE

12	NI=K1+MS-2

13	L=MIN(NI,N)

14

DO J=2,MS

15	D(J)=A(K,J)

16

ENDDO

17

DO J=K1,L

18

K2=J-K+1

19	A(K,K2)=A(K,K2)/C

20

ENDDO

21	B(K)=B(K)/C

22

23

DO I=K1,L

24

K2=I-K1+2

25

C=D(K2)

26

DO J=I,L

27

K2=J-I+1

28

K3=J-K+1

29	A(I,K2)=A(I,K2)-C*A(K,K3)

30

ENDDO

31	B(I)=B(I)-C*B(K)

32

ENDDO

33

ENDIF

34

ENDDO

35

36	IF((ABS(A(N,1))-0.000001).LE.0)THEN

37	WRITE(6,7)K

38	7     FORMAT('*****SINGULARITY IN ROW',I5)

39

STOP

40

ELSE

41	B(N)=B(N)/A(N,1)

42

DO I=1,N1

43

K=N-I

44

K1=K+1

45	NI=K1+MS-2

46	L=MIN(NI,N)

47

DO J=K1,L

48

K2=J-K+1

49	B(K)=B(K)-A(K,K2)*B(J)

50

ENDDO

51

ENDDO

52

ENDIF

53

RETURN

54

END

一个解线性方程组的程序，不知道有没有大神有空帮我看看，这是一个有限元的子程序，整个程序已经可以计算，但是结果就是不对，是一个 解二维等带宽矩阵的程序吧，数学学得不行

aliouying

麦田守望者 发表于 2014-4-18 23:04
嗯，还是谢谢你了，那样的话我还得考虑下怎么描述

如果这个程序是以前别人写的，而且是对的，那么肯定是你的矩阵组装的时候或者单元矩阵有问题
一般解方程组的这些代码是公用的，经过测试了的

如果非要检查，可以把矩阵的行号、列号、非零值 输出来，然后用matlab或者其他的求解器来对比，看求解方程组是否有问题，仅仅给定这个代码，没有算法，那是真心不好检查，只能自己猜测的事情没人愿意干。

麦田守望者

aliouying 发表于 2014-4-17 20:55
我说的输入文件的意思是指矩阵，不是其他，不要拿专业的东西，能碰到同一个专业且热心帮你解决的可能性很 ...

嗯，还是谢谢你了，那样的话我还得考虑下怎么描述

aliouying

麦田守望者 发表于 2014-4-17 19:55
恩，第一次发帖问问题，见谅哦，整体的代码有几百行，发上来也没人愿意看吧，这个应该是一个单独的方程组 ...

我说的输入文件的意思是指矩阵，不是其他，不要拿专业的东西，能碰到同一个专业且热心帮你解决的可能性很小，应该从纯数学或者fortran程序的角度上描述问题

麦田守望者

求解结果，位移，轴力，合力几乎都是错误的

麦田守望者

aliouying 发表于 2014-4-17 19:51
如果单纯给代码，咱们没法给您看，没有任何说明，没有任何算法原理

最好给出输入文件，别人才好测试你的程 ...

恩，第一次发帖问问题，见谅哦，整体的代码有几百行，发上来也没人愿意看吧，这个应该是一个单独的方程组，所以想单独拿出来，看来也不行，我把代码和输入文件发附件您帮看看行吗？

aliouying

如果单纯给代码，咱们没法给您看，没有任何说明，没有任何算法原理

最好给出输入文件，别人才好测试你的程序是否正确