107 lines
2.0 KiB
Plaintext
107 lines
2.0 KiB
Plaintext
function aitken (f, p0)
|
|
implicit none
|
|
|
|
double precision f, p0, aitken
|
|
external f
|
|
|
|
double precision p1, p2
|
|
|
|
p1 = f(p0)
|
|
p2 = f(p1)
|
|
aitken = p0 - (p1 - p0)**2 / (p2 - (2.0d0 * p1) + p0)
|
|
end
|
|
|
|
subroutine steff (f, pinit, tol, mxiter, ierr, pfinal)
|
|
implicit none
|
|
|
|
double precision f, pinit, tol, pfinal
|
|
integer mxiter, ierr
|
|
external f
|
|
|
|
double precision p0, p, aitken
|
|
integer iter
|
|
external aitken
|
|
|
|
p0 = pinit
|
|
p = aitken (f, p0)
|
|
iter = 1
|
|
while (abs (p - p0) > tol && iter < mxiter)
|
|
{
|
|
p0 = p
|
|
p = aitken (f, p0)
|
|
iter = iter + 1
|
|
}
|
|
if (abs (p - p0) > tol)
|
|
ierr = -1
|
|
else
|
|
{
|
|
ierr = 0
|
|
pfinal = p
|
|
}
|
|
end
|
|
|
|
function dcstlj (c0, c1, c2, t)
|
|
implicit none
|
|
|
|
double precision c0, c1, c2, t, dcstlj
|
|
|
|
double precision s, c01, c12, c012
|
|
|
|
s = 1.0d0 - t
|
|
c01 = (s * c0) + (t * c1)
|
|
c12 = (s * c1) + (t * c2)
|
|
c012 = (s * c01) + (t * c12)
|
|
dcstlj = c012
|
|
end
|
|
|
|
function x (t)
|
|
implicit none
|
|
double precision x, t, dcstlj
|
|
external dcstlj
|
|
x = dcstlj (2.0d0, -8.0d0, 2.0d0, t)
|
|
end
|
|
|
|
function y (t)
|
|
implicit none
|
|
double precision y, t, dcstlj
|
|
external dcstlj
|
|
y = dcstlj (1.0d0, 2.0d0, 3.0d0, t)
|
|
end
|
|
|
|
function impleq (x, y)
|
|
implicit none
|
|
double precision x, y, impleq
|
|
impleq = (5.0d0 * x * x) + y - 5.0d0
|
|
end
|
|
|
|
function f (t)
|
|
implicit none
|
|
double precision f, t, x, y, impleq
|
|
external x, y, impleq
|
|
f = impleq (x(t), y(t)) + t
|
|
end
|
|
|
|
program RCstef
|
|
implicit none
|
|
|
|
double precision x, y, impleq, f
|
|
external x, y, impleq, f, steff
|
|
|
|
double precision t0, t
|
|
integer i, ierr
|
|
|
|
t0 = 0.0d0
|
|
for (i = 0; i != 11; i = i + 1)
|
|
{
|
|
call steff (f, t0, 0.00000001d0, 1000, ierr, t)
|
|
if (ierr < 0)
|
|
write (*,*) "t0 = ", t0, " : no answer"
|
|
else if (abs (impleq (x(t), y(t))) <= 0.000001)
|
|
write (*,*) "t0 = ", t0, " : intersection at (", _
|
|
x(t), ", ", y(t), ")"
|
|
else
|
|
write (*,*) "t0 = ", t0, " : spurious solution"
|
|
t0 = t0 + 0.1d0
|
|
}
|
|
end
|