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 |
