# Program to numerically solve 3 ODE's **URL:** \ **Category:** Homework\ **Created:** [November 15, 2021, 4:54pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275 "2021-11-15T16:54:45Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 15, 2021, 4:54pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/1 "2021-11-15T16:54:46Z")
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 ![Screen Shot 2021-11-15 at 11.54.21 AM](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/c/c1735cdb3fb6d84b6ae00a19d5f641d9ce2b6c6f.png) ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [November 15, 2021, 5:13pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/2 "2021-11-15T17:13:06Z")
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. ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 15, 2021, 6:03pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/3 "2021-11-15T18:03:07Z")
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? ```auto integer, parameter :: dp = kind(1d0) ``` ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [November 15, 2021, 6:07pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/4 "2021-11-15T18:07:09Z")
No, it’s double precision 1.0. It can also be written `integer, parameter :: dp = kind(1.0d0)`, which is what I do. ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 15, 2021, 8:04pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/5 "2021-11-15T20:04:01Z")
Okay, so then I see later on in the program a variable defined later with the formatting, is this not redundant or not required? ```auto integer, parameter :: dp = kind(1d0) real(dp) :: tstart tstart = 0.0d0 ``` tstart is 0.0d0? What is going on here? RP ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [November 15, 2021, 8:42pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/6 "2021-11-15T20:42:43Z")
It is not redundant. Declaring a variable does not initialize it, so it is necessary to initialize `tstart` before using it. ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 15, 2021, 9:42pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/7 "2021-11-15T21:42:37Z")
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 ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [November 15, 2021, 10:07pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/8 "2021-11-15T22:07:23Z")
> [@RickP330](#): > > 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? 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. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [November 15, 2021, 10:18pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/9 "2021-11-15T22:18:11Z")
The d in `1.0d0` is a form of [e-notation](https://en.m.wikipedia.org/wiki/Scientific_notation#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. ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 15, 2021, 10:48pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/10 "2021-11-15T22:48:12Z")
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 ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [November 15, 2021, 11:00pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/11 "2021-11-15T23:00:45Z")
> [@RickP330](#): > > Can I have an internal function and subroutine in the same program? I’ll have both I think. Yes you can nest several subroutines and/or functions in the contains section. ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [November 15, 2021, 11:03pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/12 "2021-11-15T23:03:58Z")
> [@RickP330](#): > > 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. 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. ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [November 16, 2021, 12:54am UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/13 "2021-11-16T00:54:26Z")
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. ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 17, 2021, 7:45pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/14 "2021-11-17T19:45:22Z")
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 ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [November 17, 2021, 8:19pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/15 "2021-11-17T20:19:46Z")
One resource is the fortran-lang tutorial: [Organising code structure - Fortran Programming Language](https://fortran-lang.org/learn/quickstart/organising_code) 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](https://link.springer.com/book/10.1007/978-1-4471-6759-4). 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: ```fortran 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: ```fortran 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: ```fortran 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](https://fortran-lang.discourse.group/t/if-statement-formatting-related-to-runge-kutta-ode-solver/2240) 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. ---
**Author:** ![vsnyder](https://avatars.discourse-cdn.com/v4/letter/v/6bbea6/32.png) [@vsnyder](https://fortran-lang.discourse.group/u/vsnyder)\ **Post date:** [November 19, 2021, 11:34pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/16 "2021-11-19T23:34:38Z")
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](http://netlib.org). ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 22, 2021, 7:04pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/17 "2021-11-22T19:04:50Z")
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 ---
**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\ **Post date:** [November 22, 2021, 11:39pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/18 "2021-11-22T23:39:25Z")
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 ```auto 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 ``` ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [November 23, 2021, 8:32am UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/19 "2021-11-23T08:32:19Z")
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.) ---
**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\ **Post date:** [November 23, 2021, 3:43pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/20 "2021-11-23T15:43:47Z")
> [@RickP330](#): > > … 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. … @RickP330 , you may find it very helpful to go through [**this book by Chapman**](https://www.mheducation.com/highered/product/fortran-scientists-engineers-chapman/M9780073385891.html) first and then follow-up on online forums such as this one with targeted questions that help clarify aspects related to your particular exercise. [Next page](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275.md?page=2)