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

freedom

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

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

program mytp24
implicit none
integer i,j,n,l,k
double precision h,alpha
double precision,allocatable::a(:,:),akl(:,:),b(:),bkl(:),phi(:),x(:) ,a1(:),b1(:),c1(:),f(:),u(:)
write(6,*)'input un nombre n'
read(5,*)n
write(6,*)'input un nombre alpha'
read(5,*)alpha
allocate(a(n,n),akl(n,n),b(n),bkl(n),x(n),a1(n),b1(n),c1(n),f(n),u(n))
h=1.d0/n
x(1)=0
do i=1,n
x(i)=(i-1)*h
end do
a(1:n,1:n)=0
b(1:n)=0
b(1)=-alpha
call Integ(h,x,akl,f,n)
do i=1,n
    do j=1,n
    akl(i,j)=akl(i,j)
    bkl(i)=f(i)
enddo
end do
do l=1,2
    do k=1,n
    i=k+l-1
    j=k+l-1
if(i/=n+1.and.j/=n+1) then
a(i,j)=a(i,j)+akl(i,j)
elseif(i/=n+1) then
b(i)=b(i)+bkl(i)
endif
enddo
enddo
    do i=1,n
     write(6,'("(",i3,")=",d15.8)')i,b(i)
        do j=1,n
        write(6,'("(",i3,",",i3,")=",d15.8)')i,j,a(i,j)
        enddo
        enddo
do i=1,n
    b1(i)=a(i,i)
    enddo
    do i=2,n
    a1(i)=a(i,i-1)
    enddo
    do i=1,n-1
       c1(i)=a(i,i+1)
    end do
    do i=1,n
    f(i)=b(i)
    enddo
    a1(1)=0
    c1(n)=0
    u(n)=0
    open(25,file='le solution de u(x).txt')
    call tridag(a1,b1,c1,f,u,n)
    do i=1,n
    write(6,'("(",i3,")=",d15.8)')i,u(i)
    write(25,*)'le solution de u(x)',u(i)
    end do
    close(25)
pause
end program mytp24

subroutine Integ(h,x,akl,f,n)
implicit none
integer i,n,j
double precision h,x(n),a,b,akl(n,n),f(n),c,alpha
a=h*(1.d0/2.d0-1.d0/(2.d0*3.d0**(0.5d0)))
b=h*(1.d0/2.d0+1.d0/(2.d0*3.d0**(0.5d0)))
x(1)=0
c=1
do i=1,n
x(i)=(i-1)*h
enddo
    akl(1,1)=h/2.d0*((-(x(1)+a)/h+1)**2+((-x(1)+b)/h+1)**2)+1.d0/h
    f(1)=c*h/2.d0*((-(x(1)+a)/h+1)+((-x(1)+b)/h+1))+alpha
    do i=1,n
    f(i+1)=c*h/2.d0*((-(x(1)+a)/h+1)+((-x(1)+b)/h+1))
    do j=1,n
    if(i==j)then
    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
    elseif(i==j+1) then
    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
    elseif(i+1==j) then
    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
    else 
    akl(i,j)=0
    endif
    enddo
    enddo
    open(26,file='les matrix a et b.txt')
       do i=1,n
       write(6,'("(",i3,")=",d15.8)')i,f(i)
        do j=1,n
        write(6,'("(",i3,",",i3,")=",d15.8)')i,j,akl(i,j)
         write(26,*)'les matrix a et b',i,j,akl(i,j),i,f(i)
        enddo
        enddo 
        close(26) 
    end subroutine integ

subroutine tridag(a1,b1,c1,f,u,n)
implicit none
integer n,j
double precision gam(n),a1(n),b1(n),c1(n),u(n),f(n),bet
if (b1(1)==0.) pause 'b(1)=0 dans tridag'
 bet=b1(1)
u(1)=f(1)/bet
do j=2,n
    gam(j)=c1(j-1)/bet
    bet=b1(j)-a1(j)*gam(j)
    if (bet==0.) pause 'bet=0 dans tridag'
    u(j)=(f(j)-a1(j)*u(j-1))/bet
end do
do j=n-1,1,-1
    u(j)=u(j)-gam(j+1)*u(j+1)
end do
 end subroutine tridag

 double precision function Fbase(x,phi,dphi,n)
 integer i,j,n
 double precision phi(n,n),dphi(n,n),x(n+1)
 do i=1,n
    do j=1,n
    phi(1,1)=(-x(1)/h+1)**2
    dphi(1,1)=1.d0/h
 if(i+1==j+1)then
   phi(i,j)=(x(i)/h-i+2)**2+(x(i+1)/h+i+1)**2
   dphi(i,j)=2.d0/h
   else if(i+1==j)then
   phi(i,j)=(x(i)/h+i)*(x(j)/h+j)
   dphi(i,j)=-1.d0/h
   else if(i==j+1)then
   phi(i,j)=(x(j)/h+j)*(x(i)/h+i)
   dphi(i,j)=-1.d0/h
   else
   phi(i,j)=0
   dphi(i,j)=0
   endif
   enddo
   enddo
   end function Fbase

fcode

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

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

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

freedom

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

freedom

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