[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