Fortran : Using arrays

[SalfordPhysics]

Homework Statement

I have a code that successfully plots the trajectory of a ball moving under gravity and air resistance, but my method is rather long-winded and I want to use a 4d vector-first order ODE instead - but I don't know how to do it. I've tried writing some simple skeletons but can't get them to work either. I basically need instructions on how to do it.

Homework Equations

These are the equations that are processed to find (i+1) value of position and velocity in x and y components;
x2 = x1 + vx1*dt
y2 = y1+vy1*dt
vx2 = vx1 + Fx*dt
vy2 = vy1 + (-g*dt) + Fy*dt

And the corresponding equations for the functions;
Fx = (-1)*alpha*v*vx
Fy = (-1)*alpha*v*vy

(The x and y values are put into arrays and plotted on a graph)

My notes give Y = (x,y,vx,yv) and F(x,y,vx,vy,t) = (vx,vy,ax,ay)
With dY=dt = F(Y,t).

This part I am unsure on how to deal with.

The Attempt at a Solution

My attempts at a skeleton for the vector was as follows;
REAL, DIMENSION (1,2,3,4) :: Y

Y(1) = x1 + vx1*dt
Y(2) = y1+vy1*dt
Y(3) = vx1 + Fx*dt
Y(4) = vy1 + (-g*dt) + Fy*dt

I got rank errors at the open parentheses on all vectors. In truth, I have little idea how to do this but I know it can;t be complicated.

 

[SteamKing]

Homework Statement

I have a code that successfully plots the trajectory of a ball moving under gravity and air resistance, but my method is rather long-winded and I want to use a 4d vector-first order ODE instead - but I don't know how to do it. I've tried writing some simple skeletons but can't get them to work either. I basically need instructions on how to do it.

Homework Equations

These are the equations that are processed to find (i+1) value of position and velocity in x and y components;
x2 = x1 + vx1*dt
y2 = y1+vy1*dt
vx2 = vx1 + Fx*dt
vy2 = vy1 + (-g*dt) + Fy*dt

And the corresponding equations for the functions;
Fx = (-1)*alpha*v*vx
Fy = (-1)*alpha*v*vy

(The x and y values are put into arrays and plotted on a graph)

My notes give Y = (x,y,vx,yv) and F(x,y,vx,vy,t) = (vx,vy,ax,ay)
With dY=dt = F(Y,t).

This part I am unsure on how to deal with.

The Attempt at a Solution

My attempts at a skeleton for the vector was as follows;
REAL, DIMENSION (1,2,3,4) :: Y

Y(1) = x1 + vx1*dt
Y(2) = y1+vy1*dt
Y(3) = vx1 + Fx*dt
Y(4) = vy1 + (-g*dt) + Fy*dt

I got rank errors at the open parentheses on all vectors. In truth, I have little idea how to do this but I know it can;t be complicated.

When you declare a multi-dimensional array, each reference to an element in that array must contain the same number of dimensions as the declaration. How else can the compiler interpret which dimension or dimensions you are using if you don't specify them? Remember, you may know your intent; computers are terrible at guessing what your intent is.

 

[Mark44]

I don't think you understand what this line does:

REAL, DIMENSION (1,2,3,4) :: Y

Your thread title suggests that you intend to use multidimensional arrays, and above you are declaring Y as a four-dimensional array. However, it seems to me that all you need is a one-dimensional array with four elements, which is very different from your declaration above.

If all you need is a one-d array with four elements, use this declaration:

REAL, DIMENSION (4) :: Y

With this declaration you can access Y(1), Y(2), Y(3), and Y(4).

My notes give Y = (x,y,vx,yv) and F(x,y,vx,vy,t) = (vx,vy,ax,ay)
With dY=dt = F(Y,t).

The two references to F above don't agree. In the first line, there are four arguments, but in the second line, there are only two. It appears to me that you are confused about the different roles of arrays and functions in Fortran.

 

[SalfordPhysics]

I don't think you understand what this line does:

REAL, DIMENSION (1,2,3,4) :: Y

Your thread title suggests that you intend to use multidimensional arrays, and above you are declaring Y as a four-dimensional array. However, it seems to me that all you need is a one-dimensional array with four elements, which is very different from your declaration above.

If all you need is a one-d array with four elements, use this declaration:

REAL, DIMENSION (4) :: Y

With this declaration you can access Y(1), Y(2), Y(3), and Y(4).

Thanks Mark that's the intention on the first point precisely.

The two references to F above don't agree. In the first line, there are four arguments, but in the second line, there are only two. It appears to me that you are confused about the different roles of arrays and functions in Fortran.

There are 4 arguments + the time argument in the first reference of F. In the second, there is just Y and t, since Y has arguments vx,y,vx,vy.

I will try as you said above and come back if I can get no further.

 

[Mark44]

It would be helpful if you posted the code you have that works. That would give us some insight into what it does and some ideas on how to make it less "long-winded."

Be sure to enclose it with [ code ] and [ /code ] tags (without the spaces).

 

[SalfordPhysics]

Here it is. Not ran it since I did it in labs last week mind.

 PROGRAM eulertrajectorysimple
USE DISLIN
  IMPLICIT NONE
  INTEGER, PARAMETER :: n = 600                           !Array size
  INTEGER :: i                                            !Loop index
  REAL, DIMENSION (n) :: x, y, vx, vy, t, p, q, vp, vq    !Arrays for data
  REAL, PARAMETER :: v = 50, theta = 30                   !Initial velocity, initial angle
  REAL, PARAMETER :: pi = 3.14159

  t(1) = 0                                    !Initialization
  x(1) = 0
  y(1) = 0
  vx(1) = v*cos(pi*theta/180)              !Initial x-velocity (depends on v, theta)
  vy(1) = v*sin(pi*theta/180)              !Initial y-velocity (depends on v, theta)
  vp(1) = vx(1)
  vq(1) = vy(1)


  DO i = 1,(n-1)
      CALL evalS1(x(i),y(i),vx(i),vy(i),x(i+1),y(i+1),vx(i+1),vy(i+1),v,p(i),q(i),vp(i),vq(i), p(i+1),q(i+1),vp(i+1),vq(i+1))
      PRINT*, "x=",x(i), "y=",y(i)
      PRINT*, "vx=",vx(i), "vy=",vy(i)
  END DO
  PRINT*, vx(1), vy(1)
  CALL METAFL('XWIN')
  CALL DISINI()
  CALL GRAF(0.,100.,0.,0.5,0.,50.,0.,0.5)
  CALL COLOR ('RED')
  CALL CURVE (x,y,i)
  DO i = 1,(n-1),5
      CALL RLSYMB(3, x(i), y(i))
  END DO
  CALL COLOR('BLUE')
  CALL CURVE (p,q,i)
  DO i = 1,(n-1),5
      CALL RLSYMB(3, p(i), q(i))
  END DO

  CALL DISFIN()

END PROGRAM eulertrajectorysimple

SUBROUTINE evalS1(x1, y1, vx1, vy1, x2, y2, vx2, vy2, vinitial, p1, q1, vp1, vq1, p2,     q2, vp2, vq2)   !Drag included
      IMPLICIT NONE
      REAL :: x1, y1, vx1, vy1, vinitial    !Inputs (Drag)
      REAL :: x2, y2, vx2, vy2              !Outputs (Drag)
      REAL :: p1, q1, vp1, vq1              !Inputs (NO drag)
      REAL :: p2, q2, vp2, vq2              !Outputs (NO drag)
      REAL, PARAMETER :: dt = 0.01                
      REAL, PARAMETER :: g = 9.81               
      REAL :: Fx, Fy       !Drag force's (components)
    
      CALL evalS2(vx1, vy1, vinitial, Fx, Fy, dt)
       x2 = x1 + vx1*dt                         !Drag
       y2 = y1 + vy1*dt
       vx2 = vx1 + Fx*dt
       vy2 = vy1 + (-g*dt) + Fy*dt              !Drag
       p2 = p1 + vp1*dt                       !No drag
       q2 = q1 + vq1*dt
       vp2 = vp1
       vq2 = vq1 + ((-1)*g*dt)               !No drag
END SUBROUTINE evalS1

 

SUBROUTINE evalS2(vx, vy, v, Fx1, Fy1, dt1)         ! Drag factor
     IMPLICIT NONE
     REAL :: vx, vy, v, dt1
     REAL :: Fx1, Fy1
     REAL, PARAMETER :: alpha = 0.5                 !Drag constant

     Fx1 = (-1)*alpha*vx
     Fy1 = (-1)*alpha*vy
END SUBROUTINE evalS2

 

[Mark44]

Your evalS2 subroutine has dt1 as a formal parameter but doesn't use it. You should take it out.

Your evalS1 subroutine has 17 formal parameters. You could simplify things by passing arrays rather than individual elements of arrays. The PDF in this link -- http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/F90-Array.pdf -- talks about using arrays as parameters of subroutines on p. 30.

 

[SalfordPhysics]

Your evalS2 subroutine has dt1 as a formal parameter but doesn't use it. You should take it out.

Your evalS1 subroutine has 17 formal parameters. You could simplify things by passing arrays rather than individual elements of arrays. The PDF in this link -- http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/F90-Array.pdf -- talks about using arrays as parameters of subroutines on p. 30.

Thanks Mark I will check that out.

 

[SalfordPhysics]

The thing I am not grasping is that before I had 4 arrays, now I have one array with 4 elements, so am I right in assuming that I will only be able to overwrite the values of Y(n) on each iteration? If that's the case, then how do I go about my plot?

 

[Mark44]

The thing I am not grasping is that before I had 4 arrays, now I have one array with 4 elements, so am I right in assuming that I will only be able to overwrite the values of Y(n) on each iteration? If that's the case, then how do I go about my plot?

I'm not following you. In the code you posted in post #6, you had 9 arrays: x, y, vx, vy, t, p, q, vp, and vq. Per the declaration below, each of these arrays has n (= 600) elements

REAL, DIMENSION (n) :: x, y, vx, vy, t, p, q, vp, vq    !Arrays for data

In the x array, for example, you can access the elements as x(1), x(2), x(3), ..., x(n). The other 8 arrays are accessed in a similar manner. All of the arrays are one-dimensional inasmuch as there is one index that identifies a particular array element.

In your attempted upgrade you have this declaration:

REAL, DIMENSION (1,2,3,4) :: Y

This does NOT mean that the array Y has four elements - it has 24 elements. This is a four-dimensional array, meaning that it takes four indexes to specify a given array element. You DON'T access array elements here by Y(1), Y(2), Y(3), and Y(4). I'm almost certain that you wrote this code not understanding how arrays work.

 

[SalfordPhysics]

I'm not following you. In the code you posted in post #6, you had 9 arrays: x, y, vx, vy, t, p, q, vp, and vq. Per the declaration below, each of these arrays has n (= 600) elements

REAL, DIMENSION (n) :: x, y, vx, vy, t, p, q, vp, vq    !Arrays for data

In the x array, for example, you can access the elements as x(1), x(2), x(3), ..., x(n). The other 8 arrays are accessed in a similar manner. All of the arrays are one-dimensional inasmuch as there is one index that identifies a particular array element.

In your attempted upgrade you have this declaration:

REAL, DIMENSION (1,2,3,4) :: Y

This does NOT mean that the array Y has four elements - it has 24 elements. This is a four-dimensional array, meaning that it takes four indexes to specify a given array element. You DON'T access array elements here by Y(1), Y(2), Y(3), and Y(4). I'm almost certain that you wrote this code not understanding how arrays work.

You are totally right I wrote it with an incorrect understanding.
I said 4 arrays rather poorly assuming you would know what I meant: i had 4 arrays for a trajectory with air resistance and 4 arrays without, plus the time array.
In developing the simplification, I am first doing things only with air resistance.
It's totally correct that I was intending to create Y(4) rather than Y(1,2,3,4) - you put me right on that earlier.
This is where I am struggling...
I have Y(4), and the individual elements are Y(1), Y(2), Y(3), Y(4). Therefore, I have replaced what was previously 4 individual arrays with a single 4D vector. Therefore, I am assuming I can no longer just store up all the values within Y for plotting, and I assume I need to include some new arrays for storing up the values.
Is this correct?

 

[Mark44]

I have Y(4), and the individual elements are Y(1), Y(2), Y(3), Y(4).

Y(4) denotes the fourth element in the array. With the declaration below, you have an array Y, with four elements as listed above.

REAL, DIMENSION (4) :: Y

Therefore, I have replaced what was previously 4 individual arrays with a single 4D vector.

Here is where it gets confusing. In computer science parlance, Y is already a vector, with four elements.

A two-dimensional array would be declared like this:

REAL, DIMENSION (2, 3) :: Table

As a mathematical array it would look like this:
$$\begin{bmatrix} 3.0 & 2.0 & 0.0 \\ -1.5 & 6.2 & 3.8\end{bmatrix}$$

If my Table variable were initialized to have the same values as the matrix above, the expression Table(2, 1) would be the entry in row 2, column 1, or -1.5.

A three dimensional Fortran array could be visualized as a sequence of two-dimensional tables, and it would require three indexes to specify a single number in the three-d array: one to identify the table, one to identify the row, and one to identify the column.

Therefore, I am assuming I can no longer just store up all the values within Y for plotting, and I assume I need to include some new arrays for storing up the values.
Is this correct?

I don't think so, but it's still not clear to me what you're trying to do.

 

[SalfordPhysics]

Y(4) denotes the fourth element in the array. With the declaration below, you have an array Y, with four elements as listed above.

REAL, DIMENSION (4) :: Y

Here is where it gets confusing. In computer science parlance, Y is already a vector, with four elements.


I don't think so, but it's still not clear to me what you're trying to do.

All I am trying to do is make my program more efficient by use of Y.
I am trying to make Y(1) = x, Y(2) = y, Y(3) = vx, Y(4) = vy : where I had arrays for each of them previously.
So saying I want to plot x.vs.y, that's Y(1) vs Y(2), noting that Y(1) and Y(2) are only elements, is it the case I need more arrays to hold these variables, say, xset and yset?

 

[Mark44]

I'm not sure you need arrays at all, provided that you display the currently computed values. To get the current value of, say, x, calculate x + vx * dt, and store that in x, overwriting the old value with the new value. You don't need to store all the values of x, y, etc. if you display them each time evalS1 is called.

SUBROUTINE evalS1(x, y, vx, vy, vinitial, p, q, vp, vq)   !Drag included
      IMPLICIT NONE
      REAL :: x, y, vx, vy, vinitial    !Inputs/outputs (Drag)
      REAL :: p, q, vp, vq              !Inputs/outputs (NO drag)
      REAL, PARAMETER :: dt = 0.01               
      REAL, PARAMETER :: g = 9.81             
      REAL :: Fx, Fy       !Drag force's (components)
  
      CALL evalS2(vx, vy, vinitial, Fx, Fy, dt)
       x = x + vx * dt                         !Drag
       y = y + vy*dt
       vx = vx + Fx*dt
       vy = vy + (-g*dt) + Fy*dt              !Drag
       p = p + vp*dt                       !No drag
       q = q + vq*dt
       !vp2 = vp1
       vq = vq + ((-1)*g*dt)               !No drag
      write (*, *) ...            ! Display the values of x, y, vx, vy, etc.
END SUBROUTINE evalS1

The changes I made:
1. x1 and x2 are replaced by x. Same for y1 and y2, etc.
2. The line vp2 = vp1 is commented out. As far as I can tell, it's not needed at all, since vp2 doesn't seem to be used anywhere.
3. I added a write statement at the end of the subroutine. Since old values of variables are being overwritten, the new values of x, y, vx, vy, etc. need to be displayed each time this subroutine is called.

 

[SalfordPhysics]

I'm not sure you need arrays at all, provided that you display the currently computed values. To get the current value of, say, x, calculate x + vx * dt, and store that in x, overwriting the old value with the new value. You don't need to store all the values of x, y, etc. if you display them each time evalS1 is called.
The changes I made:
1. x1 and x2 are replaced by x. Same for y1 and y2, etc.
2. The line vp2 = vp1 is commented out. As far as I can tell, it's not needed at all, since vp2 doesn't seem to be used anywhere.
3. I added a write statement at the end of the subroutine. Since old values of variables are being overwritten, the new values of x, y, vx, vy, etc. need to be displayed each time this subroutine is called.

Thanks Mark. Regarding point 2, it is the x-velocity for the plotted trajectory that has no air resistance, so that my graph shows a comparison of the effects of air resistance.
There is only one concern I have with the points you raised - how exactly can do I go from the values being printed out as you wrote in the code to putting those values on a graph? Do I need to write them to a file?

 

[Mark44]

In the code you have in post #6, you have this loop in your main program.

DO i = 1,(n-1)
      CALL evalS1(x(i),y(i),vx(i),vy(i),x(i+1),y(i+1),vx(i+1),vy(i+1),v,p(i),q(i),vp(i),vq(i), p(i+1),q(i+1),vp(i+1),vq(i+1))
      PRINT*, "x=",x(i), "y=",y(i)
      PRINT*, "vx=",vx(i), "vy=",vy(i)
END DO

The loop above calls evalS1 599 times, and uses x(i) to calculate a value for x(i + 1), and uses y(i) to calculate a value for y(i + 1), and so on. An alternate way, and better, IMO, would be to pass arrays as parameters in evalS1 rather than individual array elements. So instead of having the do loop in the main program, the loop would be in evalS1.

The code above would look like this:

CALL evalS1(x, y, vx, vy, v, p, q, vp, vq, n)

Note that I added a parameter n, the number of elements in each array.

After evalS1 returns, the x, y, vx, vy, p, q, vp, and vq arrays would be filled, and you could do whatever plots you wanted.

The code for evalS1() would look like this:

SUBROUTINE evalS1(x, y, vx, vy, vinitial, p, q, vp, vq, n)   !Drag included
      IMPLICIT NONE
      REAL, DIMENSION(n):: x, y, vx, vy, vinitial    !Inputs/outputs (Drag)
      REAL, DIMENSION(n):: p, q, vp, vq              !Inputs/outputs (NO drag)
      REAL, PARAMETER :: dt = 0.01              
      REAL, PARAMETER :: g = 9.81            
      REAL :: Fx, Fy       !Drag force's (components)
      INTEGER n    ! Number of elements in each array
      INTEGER i

      CALL evalS2(vx, vy, vinitial, Fx, Fy, dt)
      DO i = 1, n - 1
       x(i + 1) = x(i) + vx(i) * dt                         !Drag
       y(i + 1) = y(i) + vy(i) * dt
       vx(i + 1) = vx(i) + Fx * dt
       vy(i + 1) = vy(i) + (-g * dt) + Fy * dt              !Drag
       p(i + 1) = p(i) + vp(i) * dt                       !No drag
       q(i + 1) = q(i) + vq(i)*dt
       vp(i + 1) = vp(i)
       vq(i + 1) = vq(i) + ((-1) * g  *dt)               !No drag
END SUBROUTINE evalS1

 

[SalfordPhysics]

In the code you have in post #6, you have this loop in your main program.

DO i = 1,(n-1)
      CALL evalS1(x(i),y(i),vx(i),vy(i),x(i+1),y(i+1),vx(i+1),vy(i+1),v,p(i),q(i),vp(i),vq(i), p(i+1),q(i+1),vp(i+1),vq(i+1))
      PRINT*, "x=",x(i), "y=",y(i)
      PRINT*, "vx=",vx(i), "vy=",vy(i)
END DO

The loop above calls evalS1 599 times, and uses x(i) to calculate a value for x(i + 1), and uses y(i) to calculate a value for y(i + 1), and so on. An alternate way, and better, IMO, would be to pass arrays as parameters in evalS1 rather than individual array elements. So instead of having the do loop in the main program, the loop would be in evalS1.

The code above would look like this:

CALL evalS1(x, y, vx, vy, v, p, q, vp, vq, n)

Note that I added a parameter n, the number of elements in each array.

After evalS1 returns, the x, y, vx, vy, p, q, vp, and vq arrays would be filled, and you could do whatever plots you wanted.

Thanks again Mark, almost there I think.
Keeping things simplified only considering x,y,vx and vy...

What I want to do is use

 REAL, DIMENSION (4) :: Y

for the x, y, vx and vy arrays ( so Y(1) is for x, Y(2) for y etc...
Considering that Y(1) -> Y(4) are elements, not arrays, I can't just store up the values as in your code.
Therefore the options I can see are;
(a) by continually copying the calculated values from Y(1) and Y(2) into 2 additional arrays, say xset, yset and plotting from these.
(b) writing the data into a textfile and then reading the data from that file into a graph.
(c) having arrays within an array, so Y(1) is the x array, and so on (this is somethingI haven't come across at all in fortran yet)

Could you offer more of your expert advice?

 

[Mark44]

You could declare Y like this:

REAL, DIMENSION(4, n) :: Y

The array elements Y(1, 1), Y(1, 2), ..., Y(1, n) could hold the x values.
The array elements Y(2, 1), Y(2, 2), ..., Y(2, n) could hold the y values.
The array elements Y(3, 1), Y(3, 2), ..., Y(3, n) could hold the vx values.
The array elements Y(4, 1), Y(4, 2), ..., Y(4, n) could hold the vy values.

I don't think this is a good idea, though. You have the convenience of having all your data in a single array, but this convenience comes at the cost of being harder to comprehend the meaning of a particular expression. For example, you couldn't tell what an expression Y(i, j) represented without knowing the value of i, and even then you would have to look up the difference between Y(1, *), Y(2, *) and so on. The name of the array, Y, is not suggestive of the values it holds (other than the y values in row 2). Good programming practice dictates that variable names should be chosen so that human readers can easily grasp what the variable represents. If you have an array named vx, one could guess correctly that it holds a set of x-components of velocities.

 

[SalfordPhysics]

You could declare Y like this:

REAL, DIMENSION(4, n) :: Y

The array elements Y(1, 1), Y(1, 2), ..., Y(1, n) could hold the x values.
The array elements Y(2, 1), Y(2, 2), ..., Y(2, n) could hold the y values.
The array elements Y(3, 1), Y(3, 2), ..., Y(3, n) could hold the vx values.
The array elements Y(4, 1), Y(4, 2), ..., Y(4, n) could hold the vy values.

I don't think this is a good idea, though. You have the convenience of having all your data in a single array, but this convenience comes at the cost of being harder to comprehend the meaning of a particular expression. For example, you couldn't tell what an expression Y(i, j) represented without knowing the value of i, and even then you would have to look up the difference between Y(1, *), Y(2, *) and so on. The name of the array, Y, is not suggestive of the values it holds (other than the y values in row 2). Good programming practice dictates that variable names should be chosen so that human readers can easily grasp what the variable represents. If you have an array named vx, one could guess correctly that it holds a set of x-components of velocities.

Thanks Mark.
Indeed this is more of what I was after. The purpose of the modification is for efficiency only, and since the purpose is to just plot a graph from the required data, it doesn't really matter about i and j as long as I make a note(!) alongside things for myself, that side, I don;t disagree with you at all.
And again I am in total agreement with you on the labelling, it is just the symbol used in my notes, but of course it needs changing!