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
-
1
- #3
If the numerical values of an unknown function y = f(x) are given as 2 arrays (with enough small step):
then the array of derivative could be aproximatelly computed using the formula
This is not complicated to program.
Code:
x[1],..,x[i],..,x[n]
y[1],..,y[i],..,y[n]
Code:
for i=1,..,n-1:
dy[i] = (y(i+1]-y[i])/(x(i+1]-x[i])- Thread starter
- #4
Hi, I have a data set like
x[1],..,x,..,x[n]
y[1],..,y,..,y[n]
But there is no function like y = f(x), it is unknown, and the step between the x values is not constant.
x values are calculated separately from y, but finally i have to find partial derivative like dy/dx.
Can I calculate it with
for i=1,..,n-1:
dy = (y(i+1]-y)/(x(i+1]-x)?
Or do I need to use cubic spline and then from spline parameters find the derivative? I read about it, but for me it is complicated to program such a routine.
Thank's in advance
x[1],..,x,..,x[n]
y[1],..,y,..,y[n]
But there is no function like y = f(x), it is unknown, and the step between the x values is not constant.
x values are calculated separately from y, but finally i have to find partial derivative like dy/dx.
Can I calculate it with
for i=1,..,n-1:
dy = (y(i+1]-y)/(x(i+1]-x)?
Or do I need to use cubic spline and then from spline parameters find the derivative? I read about it, but for me it is complicated to program such a routine.
Thank's in advance
modelling said:
Ok, here is a simple example, which should help you to start:
deriv.f90
Code:
[COLOR=#a020f0]program[/color] deriv
[COLOR=#0000ff]! declare variables[/color]
[COLOR=#2e8b57][b]integer[/b][/color] :: i
[COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: N [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]9[/color]
[COLOR=#2e8b57][b] real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](N) :: x, y, dy
[COLOR=#0000ff]! assign data[/color]
x [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]0[/color]., [COLOR=#ff00ff]0.1[/color], [COLOR=#ff00ff]0.2[/color], [COLOR=#ff00ff]0.3[/color], [COLOR=#ff00ff]0.5[/color], [COLOR=#ff00ff]0.6[/color], [COLOR=#ff00ff]0.8[/color], [COLOR=#ff00ff]0.9[/color], [COLOR=#ff00ff]1[/color]. [COLOR=#804040][b]/[/b][/color])
y [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]0[/color]., [COLOR=#ff00ff]0.01[/color], [COLOR=#ff00ff]0.04[/color], [COLOR=#ff00ff]0.09[/color], [COLOR=#ff00ff]0.25[/color], [COLOR=#ff00ff]0.36[/color], [COLOR=#ff00ff]0.64[/color], [COLOR=#ff00ff]0.81[/color], [COLOR=#ff00ff]1[/color]. [COLOR=#804040][b]/[/b][/color])
[COLOR=#0000ff]! compute derivation[/color]
[COLOR=#804040][b]do[/b][/color] i [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color], N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
dy(i) [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))
[COLOR=#804040][b]end do[/b][/color]
[COLOR=#0000ff]! computing last derivation using linear extrapolation[/color]
dy(N) [COLOR=#804040][b]=[/b][/color] dy(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] (dy(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]dy(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color]))[COLOR=#804040][b]/[/b][/color](x(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]x(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color]))[COLOR=#804040][b]*[/b][/color](x(N)[COLOR=#804040][b]-[/b][/color]x(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))
[COLOR=#0000ff]! print the results[/color]
[COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]10[/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]1[/color], N
[COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]20[/color]) x(i), y(i), dy(i)
[COLOR=#804040][b]end do[/b][/color]
[COLOR=#6a5acd]10[/color] [COLOR=#804040][b]format[/b][/color](a10, a10, a10)
[COLOR=#6a5acd]20[/color] [COLOR=#804040][b]format[/b][/color]([COLOR=#008080]f10.2[/color], [COLOR=#008080]f10.2[/color], [COLOR=#008080]f10.2[/color])
[COLOR=#a020f0]end program[/color] deriv
Code:
$ g95 deriv.f90 -o deriv
$ deriv
x y=f(x) dy/dx
0.00 0.00 0.10
0.10 0.01 0.30
0.20 0.04 0.50
0.30 0.09 0.80
0.50 0.25 1.10
0.60 0.36 1.40
0.80 0.64 1.70
0.90 0.81 1.90
1.00 1.00 2.10I used the function y = x*x for x from [0, 1].
You know that the derivative of this function is dy/dx = 2*x, so you can see the numerical inaccuracies.
To find the more accurate approach for computing derivation read something about numerical mathematics.
-
2
- #7
Hi all
I made a short program to check the derivative of the y=x*x function using spline. As the result showns below, the derivative (dy/dx = 2*x) is accurate except near the ends:
x y=f(x) dy/dx
.00 .00 .10
.10 .01 .17
.20 .04 .41
.30 .09 .60
.40 .16 .80
.50 .25 1.00
.60 .36 1.20
.70 .49 1.40
.80 .64 1.60
.90 .81 1.80
1.00 1.00 2.00
1.10 1.21 2.20
1.20 1.44 2.40
1.30 1.69 2.60
1.40 1.96 2.80
1.50 2.25 3.00
1.60 2.56 3.20
1.70 2.89 3.40
1.80 3.24 3.60
1.90 3.61 3.80
2.00 4.00 4.00
2.10 4.41 4.20
2.20 4.84 4.40
2.30 5.29 4.59
2.40 5.76 4.84
2.50 6.25 4.85
2.60 6.76 5.77
2.70 7.29 3.26
2.80 7.84 13.59
2.90 8.41 5.70
The program:
Program PSpline
implicit none
real*8 x(100),y(100),c(4,100)
real*8 b(100),diag(100),sub(100),sup(100)
real*8 df1,df3,dx,b1,bn
integer*4 i,n
n = 20
do i=1,n
x(i) = 0.1*(i-1)
y(i) = x(i)**2
enddo
c ----- start of spline
b(1) = (y(2)-y(1))/(x(2)-x(1))
b
= (y-y(n-1))/(x-x(n-1))
b1 = b(1)
bn = b
do i = 2,n-1
sub(i) = x(i+1)-x(i)
diag(i) = 2*(x(i+1)-x(i-1))
sup(i) = x(i)-x(i-1)
b(i) = (x(i+1)-x(i))*(y(i)-y(i-1))/(x(i)-x(i-1))
b(i) = b(i) + (x(i)-x(i-1))*(y(i+1)-y(i))/(x(i+1)-x(i))
b(i) = b(i)*3
enddo
b(2) = b(2)-(x(3)-x(2))*b(1)
b(n-1) = b(n-1)-(x(n-1)-x(n-2))*b
do i = 1,n-2
b(i) = b(i+1)
sub(i) = sub(i+1)
sup(i) = sup(i+1)
diag(i) = diag(i+1)
enddo
c ----- start trid
if ((n-2).le.1) then
b(1) = b(1)/diag(1)
goto 20
endif
do i = 2,n-2
sub(i)= sub(i)/diag(i-1)
diag(i) = diag(i)-sub(i)*sup(i-1)
b(i) = b(i)-sub(i)*b(i-1)
enddo
b(n-2) = b(n-2)/diag(n-2)
do i = n-1,1,-1
b(i) =(b(i)-sup(i)*b(i+1))/diag(i)
enddo
c ----- end trid
20 continue
do i = n-2,1,-1
b(i+1) = b(i)
enddo
b(1) = b1
do i = 1,n
c(1,i) = y(i)
c(2,i) = b(i)
enddo
do i = 1,n-1
dx = x(i+1)-x(i)
df1 = ( c(1,i+1)-c(1,i))/dx
df3 = c(2,i) + c(2,i+1)-2*df1
c(3,i) = (df1-c(2,i)-df3)/dx
c(4,i) = df3/dx**2
enddo
c ----- end of spline
write (*,'(3a10)') 'x', 'y=f(x)', 'dy/dx'
do i=1,n
write(1,'(5f10.2)') x(i), y(i), c(2,i)
enddo
999 return
end
I made a short program to check the derivative of the y=x*x function using spline. As the result showns below, the derivative (dy/dx = 2*x) is accurate except near the ends:
x y=f(x) dy/dx
.00 .00 .10
.10 .01 .17
.20 .04 .41
.30 .09 .60
.40 .16 .80
.50 .25 1.00
.60 .36 1.20
.70 .49 1.40
.80 .64 1.60
.90 .81 1.80
1.00 1.00 2.00
1.10 1.21 2.20
1.20 1.44 2.40
1.30 1.69 2.60
1.40 1.96 2.80
1.50 2.25 3.00
1.60 2.56 3.20
1.70 2.89 3.40
1.80 3.24 3.60
1.90 3.61 3.80
2.00 4.00 4.00
2.10 4.41 4.20
2.20 4.84 4.40
2.30 5.29 4.59
2.40 5.76 4.84
2.50 6.25 4.85
2.60 6.76 5.77
2.70 7.29 3.26
2.80 7.84 13.59
2.90 8.41 5.70
The program:
Program PSpline
implicit none
real*8 x(100),y(100),c(4,100)
real*8 b(100),diag(100),sub(100),sup(100)
real*8 df1,df3,dx,b1,bn
integer*4 i,n
n = 20
do i=1,n
x(i) = 0.1*(i-1)
y(i) = x(i)**2
enddo
c ----- start of spline
b(1) = (y(2)-y(1))/(x(2)-x(1))
b
= (y-y(n-1))/(x-x(n-1))b1 = b(1)
bn = b
do i = 2,n-1
sub(i) = x(i+1)-x(i)
diag(i) = 2*(x(i+1)-x(i-1))
sup(i) = x(i)-x(i-1)
b(i) = (x(i+1)-x(i))*(y(i)-y(i-1))/(x(i)-x(i-1))
b(i) = b(i) + (x(i)-x(i-1))*(y(i+1)-y(i))/(x(i+1)-x(i))
b(i) = b(i)*3
enddo
b(2) = b(2)-(x(3)-x(2))*b(1)
b(n-1) = b(n-1)-(x(n-1)-x(n-2))*b
do i = 1,n-2
b(i) = b(i+1)
sub(i) = sub(i+1)
sup(i) = sup(i+1)
diag(i) = diag(i+1)
enddo
c ----- start trid
if ((n-2).le.1) then
b(1) = b(1)/diag(1)
goto 20
endif
do i = 2,n-2
sub(i)= sub(i)/diag(i-1)
diag(i) = diag(i)-sub(i)*sup(i-1)
b(i) = b(i)-sub(i)*b(i-1)
enddo
b(n-2) = b(n-2)/diag(n-2)
do i = n-1,1,-1
b(i) =(b(i)-sup(i)*b(i+1))/diag(i)
enddo
c ----- end trid
20 continue
do i = n-2,1,-1
b(i+1) = b(i)
enddo
b(1) = b1
do i = 1,n
c(1,i) = y(i)
c(2,i) = b(i)
enddo
do i = 1,n-1
dx = x(i+1)-x(i)
df1 = ( c(1,i+1)-c(1,i))/dx
df3 = c(2,i) + c(2,i+1)-2*df1
c(3,i) = (df1-c(2,i)-df3)/dx
c(4,i) = df3/dx**2
enddo
c ----- end of spline
write (*,'(3a10)') 'x', 'y=f(x)', 'dy/dx'
do i=1,n
write(1,'(5f10.2)') x(i), y(i), c(2,i)
enddo
999 return
end
Hi gullipe,
Very interesting your program with the splines - you get the star!
Btw, a little correction:
you write the header to the screen
but the data to the file
I tried your program and got
You used equidistant data (h=0.1) computed with the known formula y = x*x
Then I tried it with non equidistant data
and I got this
The problem is that the results computed in the last 4 points are very inaccurate.
So when the function formula is unknown and you have only the data given, this approach delivers in the last 4 points very inaccurate results.
Maybe, it would be possible to create for the given data 4 additional end points (with an extrapolation) and then the splines would deliver better results.
Very interesting your program with the splines - you get the star!
Btw, a little correction:
you write the header to the screen
Code:
write ([COLOR=red]*[/color],'(3a10)') 'x', 'y=f(x)', 'dy/dx'
Code:
write([COLOR=red]1[/color],'(5f10.2)') x(i), y(i), c(2,i)I tried your program and got
Code:
x y=f(x) dy/dx
0.00 0.00 0.10
0.10 0.01 0.17
0.20 0.04 0.41
0.30 0.09 0.60
0.40 0.16 0.80
0.50 0.25 1.00
0.60 0.36 1.20
0.70 0.49 1.40
0.80 0.64 1.60
0.90 0.81 1.80
1.00 1.00 2.00
1.10 1.21 2.20
1.20 1.44 2.40
1.30 1.69 2.59
1.40 1.96 2.83
1.50 2.25 2.90
1.60 2.56 3.57
1.70 2.89 2.02
1.80 3.24 8.77
1.90 3.61 3.70You used equidistant data (h=0.1) computed with the known formula y = x*x
Then I tried it with non equidistant data
Code:
n = 9
x = (/ 0., 0.1, 0.2, 0.3, 0.5, 0.6, 0.8, 0.9, 1. /)
y = (/ 0., 0.01, 0.04, 0.09, 0.25, 0.36, 0.64, 0.81, 1. /)
Code:
x y=f(x) dy/dx
0.00 0.00 0.10
0.10 0.01 0.17
0.20 0.04 0.40
0.30 0.09 0.61
0.50 0.25 0.94
0.60 0.36 1.38
0.80 0.64 0.66
0.90 0.81 4.52
1.00 1.00 1.90The problem is that the results computed in the last 4 points are very inaccurate.
So when the function formula is unknown and you have only the data given, this approach delivers in the last 4 points very inaccurate results.
Maybe, it would be possible to create for the given data 4 additional end points (with an extrapolation) and then the splines would deliver better results.
Hi, mikrom
Thank you for the star!
In the original program, I wrote BOTH on the screen and in a file, and I have clearly been in too much hurry, when I deleted two of the four write statements (and the open statement), before I copied the program to the web site.
By the way, how do you put the code and data in these boxes with the “CODE” header and Courier font?? It looks much better that mixing text, program and output all together.
Regarding the ends: Yes, it might be possible to create extra 4-5 points with extrapolation. But it might also be possible to use some other method than cubic-spline for the last 4-5 points. For instance fit a parabola through the last points, three at a time (without any end constrains as the cubic spline requires), and calculate the slope at the center point. … Just a suggestion, don’t know if it works well.
Thank you for the star!
In the original program, I wrote BOTH on the screen and in a file, and I have clearly been in too much hurry, when I deleted two of the four write statements (and the open statement), before I copied the program to the web site.
By the way, how do you put the code and data in these boxes with the “CODE” header and Courier font?? It looks much better that mixing text, program and output all together.
Regarding the ends: Yes, it might be possible to create extra 4-5 points with extrapolation. But it might also be possible to use some other method than cubic-spline for the last 4-5 points. For instance fit a parabola through the last points, three at a time (without any end constrains as the cubic spline requires), and calculate the slope at the center point. … Just a suggestion, don’t know if it works well.
I post the code between the tagsgullipe said:
[tt][ignore]
Code:
...Read more about TGML here (or at the bottom link Process TGML)
I found the formula for three-point estimation here
and I used it. The results should be more precise than with using two-point formula.
deriv.f90
Output:
With 3-point formula the results for first and last point must be computed at other way. I used for it linear extrapolation.
and I used it. The results should be more precise than with using two-point formula.
deriv.f90
Code:
[COLOR=#a020f0]program[/color] deriv
[COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#0000ff]! declare variables[/color]
[COLOR=#2e8b57][b]integer[/b][/color] :: i
[COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: N [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]9[/color]
[COLOR=#2e8b57][b] real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](N) :: x, y, dy
[COLOR=#0000ff]! assign data[/color]
x [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]0[/color]., [COLOR=#ff00ff]0.1[/color], [COLOR=#ff00ff]0.2[/color], [COLOR=#ff00ff]0.3[/color], [COLOR=#ff00ff]0.5[/color], [COLOR=#ff00ff]0.6[/color], [COLOR=#ff00ff]0.8[/color], [COLOR=#ff00ff]0.9[/color], [COLOR=#ff00ff]1[/color]. [COLOR=#804040][b]/[/b][/color])
y [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]0[/color]., [COLOR=#ff00ff]0.01[/color], [COLOR=#ff00ff]0.04[/color], [COLOR=#ff00ff]0.09[/color], [COLOR=#ff00ff]0.25[/color], [COLOR=#ff00ff]0.36[/color], [COLOR=#ff00ff]0.64[/color], [COLOR=#ff00ff]0.81[/color], [COLOR=#ff00ff]1[/color]. [COLOR=#804040][b]/[/b][/color])
[COLOR=#0000ff]! compute derivation[/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]
dy(i) [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][COLOR=#ff00ff]1[/color])) [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=#804040][b]end do[/b][/color]
[COLOR=#0000ff]! compute first and last derivation using linear extrapolation[/color]
dy([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] dy([COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]+[/b][/color] (dy([COLOR=#ff00ff]3[/color])[COLOR=#804040][b]-[/b][/color]dy([COLOR=#ff00ff]2[/color]))[COLOR=#804040][b]/[/b][/color](x([COLOR=#ff00ff]3[/color])[COLOR=#804040][b]-[/b][/color]x([COLOR=#ff00ff]2[/color]))[COLOR=#804040][b]*[/b][/color](x([COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]x([COLOR=#ff00ff]2[/color]))
dy(N) [COLOR=#804040][b]=[/b][/color] dy(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] (dy(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]dy(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color]))[COLOR=#804040][b]/[/b][/color](x(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color]x(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color]))[COLOR=#804040][b]*[/b][/color](x(N)[COLOR=#804040][b]-[/b][/color]x(N[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))
[COLOR=#0000ff]! print the results[/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]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), dy(i)
[COLOR=#804040][b]end do[/b][/color]
[COLOR=#a020f0]end program[/color] deriv
Code:
$ g95 deriv.f90 -o deriv
$ deriv
x y=f(x) dy/dx
0.00 0.00 -0.00
0.10 0.01 0.20
0.20 0.04 0.40
0.30 0.09 0.70
0.50 0.25 0.90
0.60 0.36 1.30
0.80 0.64 1.50
0.90 0.81 1.80
1.00 1.00 2.10Hi mikrom and the others
I took this one step further and fitted every three adjadent points with a parabola. I tested the function y(x) = (1/3)*x**3 to make it a bit different from parabola:
The program is:
and the result was:
The derivative should be dy/dx = x**2, and it seems pretty accurate. Note that there are no values for the first and last points.
I took this one step further and fitted every three adjadent points with a parabola. I tested the function y(x) = (1/3)*x**3 to make it a bit different from parabola:
The program is:
Code:
Program PPolyn
implicit none
real*8 x(100),y(100),dydx(100)
real*8 dd
real*8 r1,r2,r3,r4,s1,s2
real*8 A,B,C
integer*4 i,n
n = 19
x(1) = 0.0
do i=2,n
dd = 0.1
if((i.eq.5).or.(i.eq.8)) dd = 0.2
x(i) = x(i-1)+dd
c y(i) = x(i)**2
y(i) = (1./3.)*x(i)**3
enddo
open(1,file='out.txt',status='unknown')
write (*,'(3a10)') 'x', 'y=f(x)', 'dy/dx'
write (1,'(3a10)') 'x', 'y=f(x)', 'dy/dx'
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
enddo
do i=2,n-1
write(*,'(3f10.2)') x(i), y(i), dydx(i)
write(1,'(3f10.2)') x(i), y(i), dydx(i)
enddo
endand the result was:
Code:
x y=f(x) dy/dx
.10 .00 .01
.20 .00 .04
.30 .01 .10
.50 .04 .26
.60 .07 .36
.70 .11 .50
.90 .24 .82
1.00 .33 1.00
1.10 .44 1.21
1.20 .58 1.44
1.30 .73 1.69
1.40 .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.61The derivative should be dy/dx = x**2, and it seems pretty accurate. Note that there are no values for the first and last points.
Hi all again
I also tried the program above on non-polynomial data, namely y(x) = exp(x), with the derivative y'(x) = exp(x). That gave (where "real" is the real derivative for comparison with dy/dx):
Here, I also calculated the derivative at the ends by adding two program lines into the do-loop:
I also tried the program above on non-polynomial data, namely y(x) = exp(x), with the derivative y'(x) = exp(x). That gave (where "real" is the real derivative for comparison with dy/dx):
Code:
x y=f(x) dy/dx real
.00 1.00 1.00 1.00
.10 1.11 1.11 1.11
.20 1.22 1.22 1.22
.30 1.35 1.35 1.35
.50 1.65 1.65 1.65
.60 1.82 1.83 1.82
.70 2.01 2.02 2.01
.90 2.46 2.47 2.46
1.00 2.72 2.72 2.72
1.10 3.00 3.01 3.00
1.20 3.32 3.33 3.32
1.30 3.67 3.68 3.67
1.40 4.06 4.06 4.06
1.50 4.48 4.49 4.48
1.60 4.95 4.96 4.95
1.70 5.47 5.48 5.47
1.80 6.05 6.06 6.05
1.90 6.69 6.70 6.69
2.00 7.39 7.37 7.39Here, I also calculated the derivative at the ends by adding two program lines into the do-loop:
Code:
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
enddoHi, mikrom
OK, I can explain this:
I step through all the points from 2 to n-1. For every point (i) there is one on the left (i-1) and one on the right (i+1). We can draw an explicit parabola through these three points (just as we can draw a line through two points). The equation for a parabola is y = A*x**2+B*x+C, and for each point the do-loop finds the A, B and C for the parabola that runs through the point and the two on each side. The slope at the center point (i) is y’ = 2*A*x(i) + B. Note that A, B and C are different for each step.
In the case above we are producing n-1 separate and different parabolas. But for cubic spline, on the other hand, we are producing a new function made of short pieces of third-power polynomials, which have certain constraints at the ends: same elevation and same slope, so the whole new function is differeciable two times everywhere.
A test of still another non-polynomial function, y=sin(x) with y’=cos(x) gives:
OK, I can explain this:
I step through all the points from 2 to n-1. For every point (i) there is one on the left (i-1) and one on the right (i+1). We can draw an explicit parabola through these three points (just as we can draw a line through two points). The equation for a parabola is y = A*x**2+B*x+C, and for each point the do-loop finds the A, B and C for the parabola that runs through the point and the two on each side. The slope at the center point (i) is y’ = 2*A*x(i) + B. Note that A, B and C are different for each step.
In the case above we are producing n-1 separate and different parabolas. But for cubic spline, on the other hand, we are producing a new function made of short pieces of third-power polynomials, which have certain constraints at the ends: same elevation and same slope, so the whole new function is differeciable two times everywhere.
A test of still another non-polynomial function, y=sin(x) with y’=cos(x) gives:
Code:
x y=f(x) dy/dx real
.00 .00 1.00 1.00
.10 .10 .99 1.00
.20 .20 .98 .98
.30 .30 .95 .96
.50 .48 .87 .88
.60 .56 .82 .83
.70 .64 .76 .76
.90 .78 .62 .62
1.00 .84 .54 .54
1.10 .89 .45 .45
1.20 .93 .36 .36
1.30 .96 .27 .27
1.40 .99 .17 .17
1.50 1.00 .07 .07
1.60 1.00 -.03 -.03
1.70 .99 -.13 -.13
1.80 .97 -.23 -.23
1.90 .95 -.32 -.32
2.00 .91 -.42 -.42- Thread starter
- #16
Hi, thanks for explanations.
Where can find a program where derivatives are calculated directly from a data set x and y using spline? Such approach seems to be more accurate. But my programming level is not enough to program something like that.
Thank's in advance
Where can find a program where derivatives are calculated directly from a data set x and y using spline? Such approach seems to be more accurate. But my programming level is not enough to program something like that.
Thank's in advance
Hi, modelling
I do not know where you can find that, but you are free to use the spline program above. You would have to modify that somewhat though (for input and output etc etc). But as mikrom pointed out, the spline process is very inaccurate near the ends, so I guess that you would have to switch to some other method there, for instance to the "parabola" method described just above.
I do not know where you can find that, but you are free to use the spline program above. You would have to modify that somewhat though (for input and output etc etc). But as mikrom pointed out, the spline process is very inaccurate near the ends, so I guess that you would have to switch to some other method there, for instance to the "parabola" method described just above.
Hi gullipe,
All clear !
When we have 3 known points given as x- and y-values (x1, x2, x3), (y1, y2, y3) then with the parable interpolation this system of equations must ne fullfilled:
Solving the above system we get the constants A, B, C.
All clear !
When we have 3 known points given as x- and y-values (x1, x2, x3), (y1, y2, y3) then with the parable interpolation this system of equations must ne fullfilled:
Code:
A*x1**2 + B*x1 + C = y1
A*x2**2 + B*x2 + C = y2
A*x3**2 + B*x3 + C = y3- Status
Similar threads
