1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
| program laplace2
implicit none
integer,parameter :: n=40
real*8,dimension(1:n) :: d,x
real*8,dimension(2:n) :: L
real*8,dimension(1:n-1) :: s
integer :: i
real*8 :: dx,xexact,xi,f
dx=1.d0/(n+1)
do i=1,n
d(i)=2.d0
end do
do i=2,n
L(i)=-1.d0
end do
do i=1,n-1
s(i)=-1.d0
end do
function f(y)
implicit none
p=3
pi=4.d0*atan(1.d0)
e = 1.d0
!xi=n*dx
f= (exp(-e*y)*((pi**2)*(p**2)*sin(pi*p*y)+&
2*pi*e*p*cos(pi*p*y)+e**2*sin(pi*p*y)))*dx*dx
end function f
call factolu(n,d,l,s,x(1),f,1)
call factolu(n,d,l,s,x(1),f,2)
open(1,file='resu2')
open(2,file='exacte2')
do i=1,n
write(1,*) i*dx,x(i)
write(2,*) i*dx,xexact(i*dx)
end do
close(1)
do i=1,n
write(*,*) x(i), xexact(i)
end do
end program laplace2
function xexact(x)
implicit none
integer :: p=3
real*8 :: x,xexact,pi
pi=4.d0*atan(1.d0)
xexact=exp(1.d0+x)*sin(2*p*pi*x)
xexact=x*(1.d0-x)/2
end function xexact
subroutine factolu(n,d,u,l,x,b,ityp)
integer :: n,i,ityp
real*8, dimension(1:n) :: d,x,b
real*8, dimension(2:n) ::l
real*8, dimension(1:n-1) ::u
real*8 , dimension(:),allocatable :: y
allocate(y(n))
if (ityp.eq.1) then
do i = 2,n
l(i) = l(i)/d(i-1)
d(i) = d(i)-(l(i)*u(i-1))
enddo
else if (ityp.eq.2) then
y(1) = b(1)
do i = 2,n
y(i) = b(i) - (l(i)*y(i-1))
enddo
x(n) = y(n)/d(n)
do i = n-1,1,-1
x(i) = (y(i) - u(i)*x(i+1))/d(i)
enddo
else if (ityp.eq.3) then
b(1) = x(1)*d(1) + x(2)*u(1)
b(n) = x(n-1)*l(n) + x(n)*d(n)
do i = 2,n-1
b(i) = x(i-1)*l(i) + x(i)*d(i) + x(i+1)*u(i)
enddo
endif
deallocate(y)
return
end subroutine factolu |
Partager