Tek-Tips is the largest IT community on the Internet today!
Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!
-
Congratulations Mike Lewis on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!
calculation ov derivatives from a data set 3
- Thread starter modelling
- Start date
Hi gullipe,
I tried to solve the above linear system with Vandermonde matrix manually using GEM (Gaussian Elimination method) with partial success:
I could express A in the same simple form as you but not B - my B was a little bit complicated.
But ok, I will try it again
Then I tried to use Newton's representation of interpolation polynomial
With Newton polynomial, there is not a linear system to solve, one have only to compute the divided differences
Here is the solution (I modified your program)
Output: The result are the same as yours (they should be the same, because there is only one interpolation polynomial for given points)
I tried to solve the above linear system with Vandermonde matrix manually using GEM (Gaussian Elimination method) with partial success:
I could express A in the same simple form as you but not B - my B was a little bit complicated.
But ok, I will try it again
Then I tried to use Newton's representation of interpolation polynomial
With Newton polynomial, there is not a linear system to solve, one have only to compute the divided differences
Here is the solution (I modified your program)
Code:
[COLOR=#a020f0]Program[/color] PPolyn
[COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] :: x([COLOR=#ff00ff]100[/color]),y([COLOR=#ff00ff]100[/color]),dydx([COLOR=#ff00ff]100[/color])
[COLOR=#2e8b57][b]real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] :: dd
[COLOR=#2e8b57][b]real[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]8[/color] :: fx0, fx0x1, fx1x2, fx0x1x2, a0, a1, a2, A, B
[COLOR=#2e8b57][b]integer[/b][/color][COLOR=#804040][b]*[/b][/color][COLOR=#ff00ff]4[/color] :: i,n
n [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]19[/color]
x([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0.0[/color]
[COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]2[/color],n
dd [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0.1[/color]
[COLOR=#804040][b]if[/b][/color]((i[COLOR=#804040][b].eq.[/b][/color][COLOR=#ff00ff]5[/color])[COLOR=#804040][b].or.[/b][/color](i[COLOR=#804040][b].eq.[/b][/color][COLOR=#ff00ff]8[/color])) dd [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0.2[/color]
x(i) [COLOR=#804040][b]=[/b][/color] x(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]+[/b][/color]dd
[COLOR=#0000ff]! y(i) = x(i)**2[/color]
y(i) [COLOR=#804040][b]=[/b][/color] ([COLOR=#ff00ff]1[/color].[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]3[/color].)[COLOR=#804040][b]*[/b][/color]x(i)[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]3[/color]
[COLOR=#804040][b]enddo[/b][/color]
[COLOR=#804040][b]open[/b][/color]([COLOR=#ff00ff]1[/color],[COLOR=#804040][b]file[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'out.txt'[/color],[COLOR=#804040][b]status[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]'unknown'[/color])
[COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]'(3a10)'[/color]) [COLOR=#ff00ff]'x'[/color], [COLOR=#ff00ff]'y=f(x)'[/color], [COLOR=#ff00ff]'dy/dx'[/color]
[COLOR=#804040][b]write[/b][/color] ([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]'(3a10)'[/color]) [COLOR=#ff00ff]'x'[/color], [COLOR=#ff00ff]'y=f(x)'[/color], [COLOR=#ff00ff]'dy/dx'[/color]
[COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]2[/color],n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
[COLOR=#0000ff]! divided differences[/color]
fx0 [COLOR=#804040][b]=[/b][/color]y(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])
fx0x1 [COLOR=#804040][b]=[/b][/color] (y(i)[COLOR=#804040][b]-[/b][/color]y(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))[COLOR=#804040][b]/[/b][/color](x(i)[COLOR=#804040][b]-[/b][/color]x(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))
fx1x2 [COLOR=#804040][b]=[/b][/color] (y(i[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]y(i))[COLOR=#804040][b]/[/b][/color](x(i[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]x(i))
fx0x1x2 [COLOR=#804040][b]=[/b][/color] (fx1x2 [COLOR=#804040][b]-[/b][/color] fx0x1)[COLOR=#804040][b]/[/b][/color](x(i[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]x(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))
[COLOR=#0000ff]! coefficients of Newton polynomial[/color]
[COLOR=#0000ff]! p2(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1)[/color]
a0 [COLOR=#804040][b]=[/b][/color] fx0
a1 [COLOR=#804040][b]=[/b][/color] fx0x1
a2 [COLOR=#804040][b]=[/b][/color] fx0x1x2
[COLOR=#0000ff]! coefficients of an interpolation parable y(x) = A*x*x + B*x + C[/color]
[COLOR=#0000ff]! and its derivation dy/dx = 2*A*x + B[/color]
A [COLOR=#804040][b]=[/b][/color] a2
B [COLOR=#804040][b]=[/b][/color] a1 [COLOR=#804040][b]-[/b][/color] a2[COLOR=#804040][b]*[/b][/color](x(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]+[/b][/color]x(i))
dydx(i) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color][COLOR=#804040][b]*[/b][/color]A[COLOR=#804040][b]*[/b][/color]x(i) [COLOR=#804040][b]+[/b][/color] B
[COLOR=#804040][b]if[/b][/color](i[COLOR=#804040][b].eq.[/b][/color][COLOR=#ff00ff]2[/color]) dydx(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color][COLOR=#804040][b]*[/b][/color]A[COLOR=#804040][b]*[/b][/color]x(i[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] B
[COLOR=#804040][b]if[/b][/color](i[COLOR=#804040][b].eq.[/b][/color](n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])) dydx(i[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]2[/color][COLOR=#804040][b]*[/b][/color]A[COLOR=#804040][b]*[/b][/color]x(i[COLOR=#804040][b]+[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] B
[COLOR=#804040][b]enddo[/b][/color]
[COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color],n
[COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]'(3f10.2)'[/color]) x(i), y(i), dydx(i)
[COLOR=#804040][b]write[/b][/color]([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]'(3f10.2)'[/color]) x(i), y(i), dydx(i)
[COLOR=#804040][b]enddo[/b][/color]
[COLOR=#a020f0]end[/color]
Code:
$ ppolyn_rm3
x y=f(x) dy/dx
0.00 0.00 -0.01
0.10 0.00 0.01
0.20 0.00 0.04
0.30 0.01 0.10
0.50 0.04 0.26
0.60 0.07 0.36
0.70 0.11 0.50
0.90 0.24 0.82
1.00 0.33 1.00
1.10 0.44 1.21
1.20 0.58 1.44
1.30 0.73 1.69
1.40 0.91 1.96
1.50 1.13 2.25
1.60 1.37 2.56
1.70 1.64 2.89
1.80 1.94 3.24
1.90 2.29 3.61
2.00 2.67 3.99Hi, mikrom
Very nice and a star for that. Your do loop can contain at least one less programming line than mine. Now I will take Conte & deBoor’s book “Elementary Numerical Analysis” down from the shelf and fresh up my memory of Newton formula for the interpolation polynomials, Gauss elimination algorithm and even Vandermonde matrix ...
Here is how I solved the equations by hand (for those who are interested), before I did put the solution into the program:
Very nice and a star for that. Your do loop can contain at least one less programming line than mine. Now I will take Conte & deBoor’s book “Elementary Numerical Analysis” down from the shelf and fresh up my memory of Newton formula for the interpolation polynomials, Gauss elimination algorithm and even Vandermonde matrix ...
Here is how I solved the equations by hand (for those who are interested), before I did put the solution into the program:
Code:
We have three equations:
(1) y(i-1) = A*x(i-1)**2 + B*x(i-1) + C
(2) y(i) = A*x(i)**2 + B*x(i) + C
(3) y(i+1) = A*x(i+1)**2 + B*x(i+1) + C
Let us first substract equation (2) from (1):
(4) y(i-1)-y(i) = A*(x(i-1)**2-x(i)**2) + B*(x(i-1)-x(i))
Next we substract equation (3) from (2):
(5) y(i)-y(i+1) = A*(x(i)**2-x(i+1)**2) + B*(x(i)-x(i+1))
Now we multiply equation (4) with (x(i)-x(i+1)):
(6) (y(i-1)-y(i))*(x(i)-x(i+1))
= A*(x(i-1)**2-x(i)**2)*(x(i)-x(i+1)) + B*(x(i-1)-x(i))*(x(i)-x(i+1))
... and equation (5) with (x(i-1)-x(i)):
(7) (y(i)-y(i+1))*(x(i-1)-x(i))
= A*(x(i)**2-x(i+1)**2)*(x(i-1)-x(i)) + B*(x(i)-x(i+1))*(x(i-1)-x(i))
Now we substract equation (7) from equation (6) and B disappears:
(8) (y(i-1)-y(i))*(x(i)-x(i+1)) - (y(i)-y(i+1))*(x(i-1)-x(i))
= A*[(x(i-1)**2-x(i)**2)*(x(i)-x(i+1)) - (x(i)**2-x(i+1)**2)*(x(i-1)-x(i))]
To make this more simple we put:
s1 = (x(i)**2-x(i+1)**2)*(x(i-1)-x(i))
s2 = (x(i-1)**2-x(i)**2)*(x(i)-x(i+1))
r2 = (y(i)-y(i+1))*(x(i-1)-x(i))
r1 = (y(i-1)-y(i))*(x(i)-x(i+1))
and equation (8) becomes:
r1-r2 = A*[s2-s1]
or:
A = (r1-r2)/(s2-s1)
In the same way we can find B by multiplying equation (4) with (x(i)**2-x(i+1)**2) and equation (5) with (x(i-1)**2-x(i)**2) and substract the new equations from each other to eliminate A ... etc, etc.Hi gullipe,
Thanks.
The do-loop I modified can be shortened to less lines.
But I tried to show the relationship between the divided differences (fx0, fx0x1, fx0x1x2), the coefficients of Newton's form of interpolation polynomial (a0, a1, a2) and the coefficients of general interpolation parable (A, B, C).
Thanks.
The do-loop I modified can be shortened to less lines.
But I tried to show the relationship between the divided differences (fx0, fx0x1, fx0x1x2), the coefficients of Newton's form of interpolation polynomial (a0, a1, a2) and the coefficients of general interpolation parable (A, B, C).
Hi modelling, that started this long and, to me, interesting thread.
Stop thinking about spline, and use mikrom's numerical approach, and you will end up with a very smart and short program that hopefully gives you the result that you are looking for. Although we do not know the nature of the input data you are going to use (read into the program), I am pretty sure that the result will be satisfactory. I wish you all the best in making the input part, that is not complicated to program in Fortran.
Stop thinking about spline, and use mikrom's numerical approach, and you will end up with a very smart and short program that hopefully gives you the result that you are looking for. Although we do not know the nature of the input data you are going to use (read into the program), I am pretty sure that the result will be satisfactory. I wish you all the best in making the input part, that is not complicated to program in Fortran.
Hi guys,
It is so alive, here! I have faced this problem, too, several weeks ago. My one dimensional data is irregularly distributed, so I used cubic smooth spline (in Matlab) to make an approximation of the continuous distribution, then made samples (dense enough) for derivative calculation using middle difference. The result is quite good for smooth data, but poor for discontinuous one.
I am working on two dimensional problem, right now. And I searched the previous work from different domains, the candidates are:
(1) wavelet transform;
(2) partition of unity;
... and so on
I do not have any result now.
Hope this may be helpful!
----
DONG Li
It is so alive, here! I have faced this problem, too, several weeks ago. My one dimensional data is irregularly distributed, so I used cubic smooth spline (in Matlab) to make an approximation of the continuous distribution, then made samples (dense enough) for derivative calculation using middle difference. The result is quite good for smooth data, but poor for discontinuous one.
I am working on two dimensional problem, right now. And I searched the previous work from different domains, the candidates are:
(1) wavelet transform;
(2) partition of unity;
... and so on
I do not have any result now.
Hope this may be helpful!
----
DONG Li
- Thread starter
- #28
Hi, gullipe,
Thank you for the parabola method. Could you please to help me describing the following case, see below. Here the number of 'x' values is 'n' (i=1,n), I need to calculate dy1/dx, dy2/dx and dy3/dx with parabola method. I'm confused when describing the loops for i=1,n and for 'y' values varying from j=1,3. Excuse me that I can not do such a simple case, the loop for calculating the derivatives in parabola method is i=2,nc-1 and i need 'n' values of the mentioned derivatives. And I don't know how to describe the WRITE commands for the given case.
x, y1, y2 , y3, - are given
dy1/dx dy2/dx dy3/dx - to be found and printed out
Thank you in advance.
Thank you for the parabola method. Could you please to help me describing the following case, see below. Here the number of 'x' values is 'n' (i=1,n), I need to calculate dy1/dx, dy2/dx and dy3/dx with parabola method. I'm confused when describing the loops for i=1,n and for 'y' values varying from j=1,3. Excuse me that I can not do such a simple case, the loop for calculating the derivatives in parabola method is i=2,nc-1 and i need 'n' values of the mentioned derivatives. And I don't know how to describe the WRITE commands for the given case.
x, y1, y2 , y3, - are given
dy1/dx dy2/dx dy3/dx - to be found and printed out
Thank you in advance.
Hi, modelling
Two questions first:
1) Is your data in a file? If so, what is the file name.
2) If the data is in a file, is it in this form?:
x, y11, y21, y31
x, y12, y22, y32
x, y13, y23, y33
...
x, y1n, y2n, y3n
3) Is there a comma between the columns in the file.
Two questions first:
1) Is your data in a file? If so, what is the file name.
2) If the data is in a file, is it in this form?:
x, y11, y21, y31
x, y12, y22, y32
x, y13, y23, y33
...
x, y1n, y2n, y3n
3) Is there a comma between the columns in the file.
Hi modelling,
I would create a subroutine derivation(i, x, y, dydx), where i is the index, x and y are the real input arrays and dydx is the real output array.
Then I would apply in a loop the subroutine at the all
3 arrays y1, y2, y3:
do i=2,n-1
call derivation(i, x, y1, dy1dx)
call derivation(i, x, y2, dy2dx)
call derivation(i, x, y3, dy3dx)
enddo
At the and I would write all the results into a file.
I would create a subroutine derivation(i, x, y, dydx), where i is the index, x and y are the real input arrays and dydx is the real output array.
Then I would apply in a loop the subroutine at the all
3 arrays y1, y2, y3:
do i=2,n-1
call derivation(i, x, y1, dy1dx)
call derivation(i, x, y2, dy2dx)
call derivation(i, x, y3, dy3dx)
enddo
At the and I would write all the results into a file.
Hi modelling
You actually did not answer my question: Is your input data in a file or not? ... and how is it stored there, if so? Anyway, I shall make a program for you, assuming that the data is stored in a file with comma-delimited columns. And I shall write the results in another file (in comma delimited colums also). I will also asume that y1, y2 and y3 are of equal length. Give me some days to do this.
We cannot include the "derivation" subroutines in a do loop. The do loop will be inside the subroutine, which you will call three times to fill the vectors dyXdx.
You actually did not answer my question: Is your input data in a file or not? ... and how is it stored there, if so? Anyway, I shall make a program for you, assuming that the data is stored in a file with comma-delimited columns. And I shall write the results in another file (in comma delimited colums also). I will also asume that y1, y2 and y3 are of equal length. Give me some days to do this.
We cannot include the "derivation" subroutines in a do loop. The do loop will be inside the subroutine, which you will call three times to fill the vectors dyXdx.
Hi modelling
Here is the program below. But first is a program that I made to make some data (IN.DAT) and a comparison file (COMP.DAT):
IN.DAT: This is the input file, that the program below reads:
COMP.DAT: This is file with the actual derivatives, that can be compared with OUT.DAT:
Here is the main program and the subroutine "derivation". I also put in another subroutine with mikrom's smart solution to this and with less programming lines than the other.
And finally this is the OUT.DAT file with the calculated derivatives:
Here is the program below. But first is a program that I made to make some data (IN.DAT) and a comparison file (COMP.DAT):
Code:
Program PPP
implicit none
real*8 x,y1,y2,y3,dy1,dy2,dy3
integer*4 i,n
character*1 c
c = ','
open(1,file='in.dat',status='unknown')
open(2,file='comp.dat',status='unknown')
n = 21
do i=1,n
x = (i-1)/10.
y1 = x**2
y2 = (1./3.)*x**3
y3 = dsin(x)
write (1,'(f6.3,3(a,f8.3))') x,c,y1,c,y2,c,y3
dy1 = 2*x
dy2 = x**2
dy3 = dcos(x)
write (2,'(f6.3,3(a,f8.3))') x,c,dy1,c,dy2,c,dy3
enddo
endIN.DAT: This is the input file, that the program below reads:
Code:
.000, .000, .000, .000
.100, .010, .000, .100
.200, .040, .003, .199
.300, .090, .009, .296
.400, .160, .021, .389
.500, .250, .042, .479
.600, .360, .072, .565
.700, .490, .114, .644
.800, .640, .171, .717
.900, .810, .243, .783
1.000, 1.000, .333, .841
1.100, 1.210, .444, .891
1.200, 1.440, .576, .932
1.300, 1.690, .732, .964
1.400, 1.960, .915, .985
1.500, 2.250, 1.125, .997
1.600, 2.560, 1.365, 1.000
1.700, 2.890, 1.638, .992
1.800, 3.240, 1.944, .974
1.900, 3.610, 2.286, .946
2.000, 4.000, 2.667, .909COMP.DAT: This is file with the actual derivatives, that can be compared with OUT.DAT:
Code:
.000, .000, .000, 1.000
.100, .200, .010, .995
.200, .400, .040, .980
.300, .600, .090, .955
.400, .800, .160, .921
.500, 1.000, .250, .878
.600, 1.200, .360, .825
.700, 1.400, .490, .765
.800, 1.600, .640, .697
.900, 1.800, .810, .622
1.000, 2.000, 1.000, .540
1.100, 2.200, 1.210, .454
1.200, 2.400, 1.440, .362
1.300, 2.600, 1.690, .267
1.400, 2.800, 1.960, .170
1.500, 3.000, 2.250, .071
1.600, 3.200, 2.560, -.029
1.700, 3.400, 2.890, -.129
1.800, 3.600, 3.240, -.227
1.900, 3.800, 3.610, -.323
2.000, 4.000, 4.000, -.416Here is the main program and the subroutine "derivation". I also put in another subroutine with mikrom's smart solution to this and with less programming lines than the other.
Code:
Program PPolyn
implicit none
real*8 x(100),y1(100),y2(100),y3(100)
real*8 dy1(100),dy2(100),dy3(100)
integer*4 i,n
character*1 c
c = ','
open(1,file='in.dat',status='old')
n = 1
do while(.true.)
read(1,'(4f20.0)',end=22) x(n),y1(n),y2(n),y3(n)
n = n+1
enddo
22 continue
n = n-1
close(unit=1)
! call derivation(n,x,y1,dy1)
! call derivation(n,x,y2,dy2)
! call derivation(n,x,y3,dy3)
call derivation_Mikrom(n,x,y1,dy1)
call derivation_Mikrom(n,x,y2,dy2)
call derivation_Mikrom(n,x,y3,dy3)
open(2,file='out.dat',status='unknown')
do i=1,n
write(2,'(f6.3,3(a,f8.3))') x(i),c,dy1(i),c,dy2(i),c,dy3(i)
enddo
close(unit=2)
end
! --------------------------------------------------------
Subroutine derivation(n,x,y,dydx)
implicit none
real*8 x(100),y(100),dydx(100)
real*8 r1,r2,r3,r4,s1,s2
real*8 A,B,C
integer*4 i,n
do i=2,n-1
s1 = (x(i-1)-x(i))*(x(i)**2-x(i+1)**2)
s2 = (x(i)-x(i+1))*(x(i-1)**2-x(i)**2)
r1 = (y(i-1)-y(i))*(x(i)-x(i+1))
r2 = (y(i)-y(i+1))*(x(i-1)-x(i))
r3 = (y(i-1)-y(i))*(x(i)**2-x(i+1)**2)
r4 = (y(i)-y(i+1))*(x(i-1)**2-x(i)**2)
A = (r1-r2)/(s2-s1)
B = (r3-r4)/(s1-s2)
C = y(i) - A*x(i)**2 - B*x(i)
dydx(i) = 2*A*x(i) + B
if(i.eq.2) dydx(i-1) = 2*A*x(i-1) + B
if(i.eq.(n-1)) dydx(i+1) = 2*A*x(i+1) + B
enddo
return
end
! --------------------------------------------------------
Subroutine derivation_Mikrom(n,x,y,dydx)
implicit none
real*8 x(100),y(100),dydx(100)
real*8 a1,A,B
integer*4 i,n
do i=2,n-1
a1 = (y(i)-y(i-1))/(x(i)-x(i-1))
A = ((y(i+1)-y(i))/(x(i+1)-x(i)) - a1)/(x(i+1)-x(i-1))
B = a1 - A*(x(i-1)+x(i))
dydx(i) = 2*A*x(i) + B
if(i.eq.2) dydx(i-1) = 2*A*x(i-1) + B
if(i.eq.(n-1)) dydx(i+1) = 2*A*x(i+1) + B
enddo
return
end
! --------------------------------------------------------And finally this is the OUT.DAT file with the calculated derivatives:
Code:
.000, .000, -.015, 1.005
.100, .200, .015, .995
.200, .400, .045, .980
.300, .600, .090, .950
.400, .800, .165, .915
.500, 1.000, .255, .880
.600, 1.200, .360, .825
.700, 1.400, .495, .760
.800, 1.600, .645, .695
.900, 1.800, .810, .620
1.000, 2.000, 1.005, .540
1.100, 2.200, 1.215, .455
1.200, 2.400, 1.440, .365
1.300, 2.600, 1.695, .265
1.400, 2.800, 1.965, .165
1.500, 3.000, 2.250, .075
1.600, 3.200, 2.565, -.025
1.700, 3.400, 2.895, -.130
1.800, 3.600, 3.240, -.230
1.900, 3.800, 3.615, -.325
2.000, 4.000, 4.005, -.415Hi, mikrom
I was in a too much hurry as so often before and mistook (misread) your last posting as one from "modelling". As you can see above, I put the do-loop inside the "derivation" subroutine, but after some thinking, I see that the subroutine call can be within a do-loop and then the subroutine would not contain any arrays nor a do-loop. Maybe ever smarter...? Humm.
I was in a too much hurry as so often before and mistook (misread) your last posting as one from "modelling". As you can see above, I put the do-loop inside the "derivation" subroutine, but after some thinking, I see that the subroutine call can be within a do-loop and then the subroutine would not contain any arrays nor a do-loop. Maybe ever smarter...? Humm.
- Thread starter
- #34
- Thread starter
- #36
Hi gullipe and mikrom,
Yes I use an inpit file in a form x, y1,y2,y3 what you have created.
!!!!!!!!!!!!!!!!!! Input file:
.000, .000, .000, .000
.100, .010, .000, .100
.200, .040, .003, .199
.300, .090, .009, .296
.400, .160, .021, .389
.500, .250, .042, .479
.600, .360, .072, .565
.700, .490, .114, .644
.800, .640, .171, .717
.900, .810, .243, .783
1.000, 1.000, .333, .841
1.100, 1.210, .444, .891
1.200, 1.440, .576, .932
1.300, 1.690, .732, .964
1.400, 1.960, .915, .985
1.500, 2.250, 1.125, .997
1.600, 2.560, 1.365, 1.000
1.700, 2.890, 1.638, .992
1.800, 3.240, 1.944, .974
1.900, 3.610, 2.286, .946
2.000, 4.000, 2.667, .909
!!!!!!!!!!!!!!!! Then your program:
Program PPolyn
implicit none
real*8 x(100),y1(100),y2(100),y3(100)
real*8 dy1(100),dy2(100),dy3(100)
integer*4 i,n
character*1 c
c = ','
open(1,file='in.dat',status='old')
n = 1
do while(.true.)
read(1,'(4f20.0)',end=22) x,y1,y2,y3
n = n+1
enddo
22 continue
n = n-1
close(unit=1)
! call derivation(n,x,y1,dy1)
! call derivation(n,x,y2,dy2)
! call derivation(n,x,y3,dy3)
call derivation_Mikrom(n,x,y1,dy1)
call derivation_Mikrom(n,x,y2,dy2)
call derivation_Mikrom(n,x,y3,dy3)
open(2,file='out.dat',status='unknown')
do i=1,n
write(2,'(f6.3,3(a,f8.3))') x(i),c,dy1(i),c,dy2(i),c,dy3(i)
enddo
close(unit=2)
end
! --------------------------------------------------------
Subroutine derivation(n,x,y,dydx)
implicit none
real*8 x(100),y(100),dydx(100)
real*8 r1,r2,r3,r4,s1,s2
real*8 A,B,C
integer*4 i,n
do i=2,n-1
s1 = (x(i-1)-x(i))*(x(i)**2-x(i+1)**2)
s2 = (x(i)-x(i+1))*(x(i-1)**2-x(i)**2)
r1 = (y(i-1)-y(i))*(x(i)-x(i+1))
r2 = (y(i)-y(i+1))*(x(i-1)-x(i))
r3 = (y(i-1)-y(i))*(x(i)**2-x(i+1)**2)
r4 = (y(i)-y(i+1))*(x(i-1)**2-x(i)**2)
A = (r1-r2)/(s2-s1)
B = (r3-r4)/(s1-s2)
C = y(i) - A*x(i)**2 - B*x(i)
dydx(i) = 2*A*x(i) + B
if(i.eq.2) dydx(i-1) = 2*A*x(i-1) + B
if(i.eq.(n-1)) dydx(i+1) = 2*A*x(i+1) + B
enddo
return
end
! --------------------------------------------------------
Subroutine derivation_Mikrom(n,x,y,dydx)
implicit none
real*8 x(100),y(100),dydx(100)
real*8 a1,A,B
integer*4 i,n
do i=2,n-1
a1 = (y(i)-y(i-1))/(x(i)-x(i-1))
A = ((y(i+1)-y(i))/(x(i+1)-x(i)) - a1)/(x(i+1)-x(i-1))
B = a1 - A*(x(i-1)+x(i))
dydx(i) = 2*A*x(i) + B
if(i.eq.2) dydx(i-1) = 2*A*x(i-1) + B
if(i.eq.(n-1)) dydx(i+1) = 2*A*x(i+1) + B
enddo
return
end
!!!!!!!!!!!!!! And the output file:
0.000, 0.000, -0.015, 1.005
0.100, 0.200, 0.015, 0.995
0.200, 0.400, 0.045, 0.980
0.300, 0.600, 0.090, 0.950
0.400, 0.800, 0.165, 0.915
0.500, 1.000, 0.255, 0.880
0.600, 1.200, 0.360, 0.825
0.700, 1.400, 0.495, 0.760
0.800, 1.600, 0.645, 0.695
0.900, 1.800, 0.810, 0.620
1.000, 2.000, 1.005, 0.540
1.100, 2.200, 1.215, 0.455
1.200, 2.400, 1.440, 0.365
1.300, 2.600, 1.695, 0.265
1.400, 2.800, 1.965, 0.165
1.500, 3.000, 2.250, 0.075
1.600, 3.200, 2.565, -0.025
1.700, 3.400, 2.895, -0.130
1.800, 3.600, 3.240, -0.230
1.900, 3.800, 3.615, -0.325
2.000, 4.000, 4.005, -0.415
I will try what mikrom proposed with. use a loop to call subroutine:
do i=2,n-1
call derivation(i, x, y1, dy1dx)
call derivation(i, x, y2, dy2dx)
call derivation(i, x, y3, dy3dx)
enddo
Thanks a lot, you are helping me so much.
Yes I use an inpit file in a form x, y1,y2,y3 what you have created.
!!!!!!!!!!!!!!!!!! Input file:
.000, .000, .000, .000
.100, .010, .000, .100
.200, .040, .003, .199
.300, .090, .009, .296
.400, .160, .021, .389
.500, .250, .042, .479
.600, .360, .072, .565
.700, .490, .114, .644
.800, .640, .171, .717
.900, .810, .243, .783
1.000, 1.000, .333, .841
1.100, 1.210, .444, .891
1.200, 1.440, .576, .932
1.300, 1.690, .732, .964
1.400, 1.960, .915, .985
1.500, 2.250, 1.125, .997
1.600, 2.560, 1.365, 1.000
1.700, 2.890, 1.638, .992
1.800, 3.240, 1.944, .974
1.900, 3.610, 2.286, .946
2.000, 4.000, 2.667, .909
!!!!!!!!!!!!!!!! Then your program:
Program PPolyn
implicit none
real*8 x(100),y1(100),y2(100),y3(100)
real*8 dy1(100),dy2(100),dy3(100)
integer*4 i,n
character*1 c
c = ','
open(1,file='in.dat',status='old')
n = 1
do while(.true.)
read(1,'(4f20.0)',end=22) x
,y1,y2,y3n = n+1
enddo
22 continue
n = n-1
close(unit=1)
! call derivation(n,x,y1,dy1)
! call derivation(n,x,y2,dy2)
! call derivation(n,x,y3,dy3)
call derivation_Mikrom(n,x,y1,dy1)
call derivation_Mikrom(n,x,y2,dy2)
call derivation_Mikrom(n,x,y3,dy3)
open(2,file='out.dat',status='unknown')
do i=1,n
write(2,'(f6.3,3(a,f8.3))') x(i),c,dy1(i),c,dy2(i),c,dy3(i)
enddo
close(unit=2)
end
! --------------------------------------------------------
Subroutine derivation(n,x,y,dydx)
implicit none
real*8 x(100),y(100),dydx(100)
real*8 r1,r2,r3,r4,s1,s2
real*8 A,B,C
integer*4 i,n
do i=2,n-1
s1 = (x(i-1)-x(i))*(x(i)**2-x(i+1)**2)
s2 = (x(i)-x(i+1))*(x(i-1)**2-x(i)**2)
r1 = (y(i-1)-y(i))*(x(i)-x(i+1))
r2 = (y(i)-y(i+1))*(x(i-1)-x(i))
r3 = (y(i-1)-y(i))*(x(i)**2-x(i+1)**2)
r4 = (y(i)-y(i+1))*(x(i-1)**2-x(i)**2)
A = (r1-r2)/(s2-s1)
B = (r3-r4)/(s1-s2)
C = y(i) - A*x(i)**2 - B*x(i)
dydx(i) = 2*A*x(i) + B
if(i.eq.2) dydx(i-1) = 2*A*x(i-1) + B
if(i.eq.(n-1)) dydx(i+1) = 2*A*x(i+1) + B
enddo
return
end
! --------------------------------------------------------
Subroutine derivation_Mikrom(n,x,y,dydx)
implicit none
real*8 x(100),y(100),dydx(100)
real*8 a1,A,B
integer*4 i,n
do i=2,n-1
a1 = (y(i)-y(i-1))/(x(i)-x(i-1))
A = ((y(i+1)-y(i))/(x(i+1)-x(i)) - a1)/(x(i+1)-x(i-1))
B = a1 - A*(x(i-1)+x(i))
dydx(i) = 2*A*x(i) + B
if(i.eq.2) dydx(i-1) = 2*A*x(i-1) + B
if(i.eq.(n-1)) dydx(i+1) = 2*A*x(i+1) + B
enddo
return
end
!!!!!!!!!!!!!! And the output file:
0.000, 0.000, -0.015, 1.005
0.100, 0.200, 0.015, 0.995
0.200, 0.400, 0.045, 0.980
0.300, 0.600, 0.090, 0.950
0.400, 0.800, 0.165, 0.915
0.500, 1.000, 0.255, 0.880
0.600, 1.200, 0.360, 0.825
0.700, 1.400, 0.495, 0.760
0.800, 1.600, 0.645, 0.695
0.900, 1.800, 0.810, 0.620
1.000, 2.000, 1.005, 0.540
1.100, 2.200, 1.215, 0.455
1.200, 2.400, 1.440, 0.365
1.300, 2.600, 1.695, 0.265
1.400, 2.800, 1.965, 0.165
1.500, 3.000, 2.250, 0.075
1.600, 3.200, 2.565, -0.025
1.700, 3.400, 2.895, -0.130
1.800, 3.600, 3.240, -0.230
1.900, 3.800, 3.615, -0.325
2.000, 4.000, 4.005, -0.415
I will try what mikrom proposed with. use a loop to call subroutine:
do i=2,n-1
call derivation(i, x, y1, dy1dx)
call derivation(i, x, y2, dy2dx)
call derivation(i, x, y3, dy3dx)
enddo
Thanks a lot, you are helping me so much.
- Status
Similar threads
