# Program to numerically solve 3 ODE's

Hi all,
I am tasked with writing a program to solve Stagnation point flow (Hiemenz) flow. The formulation of the solution requires us to guess at a value for f’’ and watch for boundary conditions to converge.

I’ve already received so much help here and have enough snips of code running that I’ve proved to myself it can be done - but now I need to write my own code and feel ready to do it. If you guys don’t mind and are interested I’d like to post my progress here and if you can comment I’d be most thankful.

I haven’t really programmed since Fortran 77 and some of the variable format codes look intimidating. So I’d like to ask for help along the way.

As part of the solution I have chosen to use a “Bisection” method to hone in on the correct value of f’’ and that must be part of the program. To get my thoughts in order I wrote a flow chart, for those curious I’l post it below.

Much Thanks in advance and thank you for following.
Kind regards,
Rick

2 Likes

Go ahead and post your progress here. Yes, the “shooting method” (using bisection to converge the BC) is a common method to solve boundary value problems.

Hey Thanks a lot - Shooting method = Bisection okay!

Some Code I looked at had the following format descriptor. I’m kind of scratching my head here. Isn’t “1d0” just an integer?

``````integer, parameter :: dp = kind(1d0)
``````
1 Like

No, it’s double precision 1.0. It can also be written `integer, parameter :: dp = kind(1.0d0)`, which is what I do.

2 Likes

Okay, so then I see later on in the program a variable defined later with the formatting, is this not redundant or not required?

``````integer, parameter :: dp = kind(1d0)
real(dp) :: tstart
tstart = 0.0d0
``````

tstart is 0.0d0? What is going on here?
RP

It is not redundant. Declaring a variable does not initialize it, so it is necessary to initialize `tstart` before using it.

I am really struggling with the syntax of the formatting. 1d0 says to use one digit with zero digits to the right of the decimal? then what does tstart = 0.0d0 mean?
RP

No, `1d0` is just an example of a double precision constant. Any double precision constant could have used as the argument of `kind`. Note that

`print*,kind(1d0),kind(1.0d0),kind(1.000d0)`

will print the same integer (8 for gfortran and Intel Fortran) thrice. Similarly,

`print*,kind(1e0),kind(1.0e0),kind(1.000e0)`

will print the `kind` corresponding to single precision (4 for gfortran and Intel Fortran) thrice.

1 Like

The d in `1.0d0` is a form of e-notation used to represent “times ten raised to the power of”.

For single precision real literal constants, the letter used is e, while for double precision reals, d is used instead.

1 Like

Oh, excellent, Thank you Beliavsky and Ivanpribec!

Next topic, subroutine structure, This is a rather simple program, I would like to avoid having to call any external subroutines.

Can I have an internal function and subroutine in the same program? I’ll have both I think.

The Contains function seems really neat, but it also appears a little advanced. I’d like to write something on the simple side.

What structure do you like?
RP

Yes you can nest several subroutines and/or functions in the contains section.

1 Like

Yes, you can, but internal functions and subroutines inherit all the variables defined in the main program. To increase modularity it is generally better to define procedures in a separate module and `use` the module in the main program to get access to those procedures.

2 Likes

The test “g = 0” in your flow chart is too strict for a finite-precision calculation. Replace it with something reasonable such as |g| < 10-8.

1 Like

Hi Ivan,
I am struggling to understand the “Contains” section. Unlike a Subroutine or a Function which I have to interface in the beginning of a program, a Contains section I do not?
Can you recommend a good reference for program branching so I can study the options?
Much Thanks
Rick

One resource is the fortran-lang tutorial: Organising code structure - Fortran Programming Language

If available to you, I would recommend the latest edition of “Modern Fortran explained” by Metcalf, Reid & Cohen, chapter 5 - Program units and procedures.

A second option would be “Guide to Fortran 2008 Programming” by Walter Brainerd: Guide to Fortran 2008 Programming | SpringerLink. Perhaps you can access the digital version with your university account.

The keyword `contains` is simply used to collect internal subprograms (subroutines and functions) in a main program or module:

``````program main
implicit none

! ... declarations ...

! ... main program unit commands ...

contains

! ... internal subprograms ...

subroutine runge(...)

end subroutine

function dydt(t,y)

end function

end program
``````

Instead of placing the procedures in the main program, you can place them in a module:

``````module flow_solver_routines
implicit none

! ... declarations ...

contains

! ... internal subprograms ...

subroutine runge(...)

end subroutine

function dydt(t,y)

end function

end module
``````

that you then use in the main program as follows:

``````program main
use flow_solver_routines, only: runge, dydt
implicit none

! ...

call runge(dydt,tstart,tend,y)

end program
``````

In fact even subroutines and functions can have procedures nested in their own `contains` section.

On the other hand, your old program used so-called `external` subprograms, i.e. subprograms that are not located within a program unit or module. External procedures were the norm in F77, but since F90 modules are the preferred approach to organize your (sub)programs.

When using internal subprograms the compiler can check you are calling the function or subroutine correctly since the interface is visible. Internal subprograms also have access to the objects in their host program, which can be useful for many things.

2 Likes

The literature of existing software to solve boundary-value problems is extensive. If this is a homework problem, then I guess you need to attack the problem yourself. Shooting is one way. Multiple shooting, where the code gives up and restarts when errors look too large, combined with a least-squares method to minimize the “jumps” at restarts, and collocation, are other methods. If your problem is a Sturm-Liouville problem, and you need its eigenvalues, I recommend sleign (or dleign).

But if it’s a “real” problem you should just get an existing code. Start at netlib.org.

1 Like

vsnyder,
This is a programming exercise, so I’ll have to write something in my own code. But thank you for the advice, you are right on.
Rick

1 Like

Hi Ivan,
Okay, Can I do the following, it doesn’t seem to work. I want to put the use the function within the subroutine. Is there something fundamentally wrong with this?
Rick

``````program main
implicit none

! ... declarations ...

! ... main program unit commands ...

Call runge

contains

subroutine runge(...)

dydt(t,y)

end subroutine

function dydt(t,y)

end function

end program
``````

Assuming the calling sequences are correct, it looks okay to me. Are you passing the function dydt as a dummy variable to the subroutine, or is it just found directly since it’s located in the same program scope? (The function call should probably be assigned to a left-hand side.)

@RickP330 , you may find it very helpful to go through this book by Chapman first and then follow-up on online forums such as this one with targeted questions that help clarify aspects related to your particular exercise.

1 Like