老师可不可以棒忙看一下，不知道哪里的原因但是结果总...

freedom

这是用有限元法计算一维常系数常微分方程，不知道哪里出错，但是结果跟我用手算的不一样，谢谢老师

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

查看代码

打印?

001	program mytp24

002	implicit none

003	integer i,j,n,l,k

004	double precision h,alpha

005	double precision,allocatable::a(:,:),akl(:,:),b(:),bkl(:),phi(:),x(:) ,a1(:),b1(:),c1(:),f(:),u(:)

006	write(6,*)'input un nombre n'

007	read(5,*)n

008	write(6,*)'input un nombre alpha'

009	read(5,*)alpha

010	allocate(a(n,n),akl(n,n),b(n),bkl(n),x(n),a1(n),b1(n),c1(n),f(n),u(n))

011

h=1.d0/n

012

x(1)=0

013

do i=1,n

014	x(i)=(i-1)*h

015

end do

016	a(1:n,1:n)=0

017

b(1:n)=0

018	b(1)=-alpha

019	call Integ(h,x,akl,f,n)

020

do i=1,n

021

do j=1,n

022	akl(i,j)=akl(i,j)

023	bkl(i)=f(i)

024

enddo

025

end do

026

do l=1,2

027

do k=1,n

028

i=k+l-1

029

j=k+l-1

030	if(i/=n+1.and.j/=n+1) then

031	a(i,j)=a(i,j)+akl(i,j)

032	elseif(i/=n+1) then

033	b(i)=b(i)+bkl(i)

034

endif

035

enddo

036

enddo

037

do i=1,n

038	write(6,'("(",i3,")=",d15.8)')i,b(i)

039

do j=1,n

040	write(6,'("(",i3,",",i3,")=",d15.8)')i,j,a(i,j)

041

enddo

042

enddo

043

do i=1,n

044	b1(i)=a(i,i)

045

enddo

046

do i=2,n

047	a1(i)=a(i,i-1)

048

enddo

049	do i=1,n-1

050	c1(i)=a(i,i+1)

051

end do

052

do i=1,n

053

f(i)=b(i)

054

enddo

055

a1(1)=0

056

c1(n)=0

057

u(n)=0

058	open(25,file='le solution de u(x).txt')

059	call tridag(a1,b1,c1,f,u,n)

060

do i=1,n

061	write(6,'("(",i3,")=",d15.8)')i,u(i)

062	write(25,*)'le solution de u(x)',u(i)

063

end do

064

close(25)

065

pause

066	end program mytp24

067

068	subroutine Integ(h,x,akl,f,n)

069	implicit none

070	integer i,n,j

071	double precision h,x(n),a,b,akl(n,n),f(n),c,alpha

072	a=h*(1.d0/2.d0-1.d0/(2.d0*3.d0**(0.5d0)))

073	b=h*(1.d0/2.d0+1.d0/(2.d0*3.d0**(0.5d0)))

074

x(1)=0

075

c=1

076

do i=1,n

077	x(i)=(i-1)*h

078

enddo

079	akl(1,1)=h/2.d0*((-(x(1)+a)/h+1)**2+((-x(1)+b)/h+1)**2)+1.d0/h

080	f(1)=c*h/2.d0*((-(x(1)+a)/h+1)+((-x(1)+b)/h+1))+alpha

081

do i=1,n

082	f(i+1)=c*h/2.d0*((-(x(1)+a)/h+1)+((-x(1)+b)/h+1))

083

do j=1,n

084	if(i==j)then

085	akl(i,j)=h/2.d0*(((x(i-1)+a)/h-i+2)**2+((x(i-1)+b)/h-i+2)**2)+h/2.d0*((-(x(i)+a)/h+i)**2+((-x(i)+b)/h+i)**2)+2.d0/h

086	elseif(i==j+1) then

087	akl(i,j)=h/2.d0*(-(x(j)+a)/h+i)*(-(x(j)+a)/h+j)+h/2.d0*(-(x(j)+b)/h+i)*(-(x(j)+b)/h+j)-1.d0/h

088	elseif(i+1==j) then

089	akl(i,j)=h/2.d0*(-(x(i)+a)/h+i)*(-(x(i)+a)/h+j)+h/2.d0*(-(x(i)+b)/h+i)*(-(x(i)+b)/h+j)-1.d0/h

090

else

091	akl(i,j)=0

092

endif

093

enddo

094

enddo

095	open(26,file='les matrix a et b.txt')

096

do i=1,n

097	write(6,'("(",i3,")=",d15.8)')i,f(i)

098

do j=1,n

099	write(6,'("(",i3,",",i3,")=",d15.8)')i,j,akl(i,j)

100	write(26,*)'les matrix a et b',i,j,akl(i,j),i,f(i)

101

enddo

102

enddo

103

close(26)

104	end subroutine integ

105

106	subroutine tridag(a1,b1,c1,f,u,n)

107	implicit none

108	integer n,j

109	double precision gam(n),a1(n),b1(n),c1(n),u(n),f(n),bet

110	if (b1(1)==0.) pause 'b(1)=0 dans tridag'

111

bet=b1(1)

112	u(1)=f(1)/bet

113

do j=2,n

114	gam(j)=c1(j-1)/bet

115	bet=b1(j)-a1(j)*gam(j)

116	if (bet==0.) pause 'bet=0 dans tridag'

117	u(j)=(f(j)-a1(j)*u(j-1))/bet

118

end do

119	do j=n-1,1,-1

120	u(j)=u(j)-gam(j+1)*u(j+1)

121

end do

122	end subroutine tridag

123

124	double precision function Fbase(x,phi,dphi,n)

125	integer i,j,n

126	double precision phi(n,n),dphi(n,n),x(n+1)

127

do i=1,n

128

do j=1,n

129	phi(1,1)=(-x(1)/h+1)**2

130	dphi(1,1)=1.d0/h

131	if(i+1==j+1)then

132	phi(i,j)=(x(i)/h-i+2)**2+(x(i+1)/h+i+1)**2

133	dphi(i,j)=2.d0/h

134	else if(i+1==j)then

135	phi(i,j)=(x(i)/h+i)*(x(j)/h+j)

136	dphi(i,j)=-1.d0/h

137	else if(i==j+1)then

138	phi(i,j)=(x(j)/h+j)*(x(i)/h+i)

139	dphi(i,j)=-1.d0/h

140

else

141	phi(i,j)=0

142	dphi(i,j)=0

143

endif

144

enddo

145

enddo

146	end function Fbase

fcode

计算结果不一致的问题，别人帮不了你。
因为，别人不知道什么才是你希望的，正确的结果。

对于编译器和代码来说，只要能算完，都认为是正确的。不符合你的预期，那是另外的问题。
我们并不知道你的预期，所以无法帮你。

不过，你可以自行学习如何 debug 单步调试：详见：
http://debug.w.fcode.cn
http://v.fcode.cn/video-debugger.html

freedom

我昨天弄好了，少加了一个循环，谢谢老师，单步调试听起来很酷的样子，马上看起来，谢谢老师每次都回的这么及时

freedom

单步debug 好东西啊，以前只看了课程可能没看仔细这些，（因为我们老师只丢了好多材料不知道从哪下手，所以就看了你们的视频）感恩