# Program to numerically solve 3 ODE's

**URL:** <https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275>\
**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:** 15\
**Page:** 2

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [November 23, 2021, 10:20pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/21 "2021-11-23T22:20:10Z")

</div>

Okay, I beat my head upon a wall for a while on this one. Why does this work:

```auto
y = y+(k1+2*(k2+k3)+k4)/6

```

But this does not?

```auto
y+(1/6)*(k1+(2*k2)+(2*k3)+k4)

```

---

<div class="post-metadata">

**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 23, 2021, 10:43pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/22 "2021-11-23T22:43:48Z")

</div>

> [@RickP330](#):
>
> `y+(1/6)*(k1+(2*k2)+(2*k3)+k4)`

Because you have to write it like this:

```auto
y+(1._dp/6)*(k1+(2*k2)+(2*k3)+k4)

```

Where `dp` is the kind of your (double precision?) number. If it is just single precision, you can write it as:

```auto
y+(1./6)*(k1+(2*k2)+(2*k3)+k4)

```

If you write just `1/6`, the answer is 0.

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [November 24, 2021, 11:51am UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/23 "2021-11-24T11:51:49Z")

</div>

It is important to use informative compiler options, especially when debugging. For the program

```auto
integer :: k
k = 3
print*,(1/6)*k
end

```

`gfortran -Wall -Wextra xdiv.f90` says

```auto
xdiv.f90:3:9:

 print*,(1/6)*k
         1
Warning: Integer division truncated to constant '0' at (1) [-Winteger-division]

```

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [November 24, 2021, 9:46pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/24 "2021-11-24T21:46:10Z")

</div>

Got it, thank you guys…  
RP

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [November 24, 2021, 11:13pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/25 "2021-11-24T23:13:32Z")

</div>

Hi Guys,  
Here is my first working pass on the program. I think I should add some code to prevent an infinite loop if you input bad values for the guess, but if you put them in right (1.3 and 1.1 for instance). The program works.  
If you can be so kind as to point out ways to tighten up the code I’d really appreciate it.

```auto
Program AE6009
	implicit none
!
	integer, parameter :: dp = kind(1d0)
	integer, parameter :: n = 3 !Order of ODE's
	real(dp) :: a, b, c, g !Variables for shooting method
	real(dp) :: eta, deta, etastart, etastop !Variables for Runge-cutta
	real(dp), dimension(n) :: y, k1, k2, k3, k4
!	
	Write (6,*)''
	write (6,*) 'Enter a guess for f"(0) < 0'
	Read (*,*) a
	Write (6,*) 'Enter a guess for f"(0) > 0'
	Read (*,*) b
!	
	g = 1.0d0
	do while (g > 1d-10 .or. g < -1d-10)
!
	C=0.5*(a+b)
	etastart = 0.0d0
	etastop = 4.8d0 ! eta at infinity
	deta = 0.01d0
	y(1) = 0.0d0 ! Initial Condition f(0) = 0
	y(2) = 0.0d0 ! Initial Condition f'(0) = 0
	y(3) = c ! Initial Condition guess for f''(0)
	eta = etastart
!
		do while(eta < etastop)
			call Runge ()
			if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
			end if
		end do
!
	g = y(2)-1d0
!
		if (g > 0d0)then
			a = c
			else
			b = c
		end if
	end do
	write (6,*)''
	Write (6,*) 'The value for f"(0) that satisfies fprime(inf) = 1'
	Write (6,'(a,f10.5)') 'f"(0) =' ,c
	Write (6,*)''
!
! Now calculate the solutions to the ODE's
!
	etastart = 0.0d0
	etastop = 4.8d0 ! eta at infinity
	deta = 0.01d0
	y(1) = 0.0d0 ! Initial Condition f(0) = 0
	y(2) = 0.0d0 ! Initial Condition f'(0) = 0
	y(3) = c ! Initial Condition f"(0) = c
	eta = etastart
	Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
	do while(eta < etastop)
			call Runge ()
		if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
			Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
		end if
	end do
!
!	Program subroutines and sub functions 
!
Contains
	Subroutine Runge ()
		k1 = deta*f(eta,y) !Runge-cutta 4th order formulas
		k2 = deta*f(eta+deta/2, y+k1/2)
		k3 = deta*f(eta+deta/2, y+k2/2)
		k4 = deta*f(eta+deta, y+k3)
		y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
		eta = eta+deta
	end subroutine
!
	Function f(eta,y)
		real(dp), intent(in) :: eta, y(n)
		real (dp) :: f(n)
		f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
	end Function
!
end Program AE6009

```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [November 25, 2021, 1:12am UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/26 "2021-11-25T01:12:46Z")

</div>

> [@RickP330](#):
>
> If you can be so kind as to point out ways to tighten up the code I’d really appreciate it.

It looks good, but I have inserted some suggestions and questions as comments prefaced by `!!`.

```auto
Program AE6009
	implicit none
!
	integer, parameter :: dp = kind(1d0)
	integer, parameter :: n = 3 !Order of ODE's
	real(dp) :: a, b, c, g !Variables for shooting method
	real(dp) :: eta, deta, etastart, etastop !Variables for Runge-cutta
	real(dp), dimension(n) :: y, k1, k2, k3, k4
!	
	Write (6,*)'' !! use * instead of 6 or output_unit from iso_fortran_env
	write (6,*) 'Enter a guess for f"(0) < 0'
	Read (*,*) a
	Write (6,*) 'Enter a guess for f"(0) > 0'
	Read (*,*) b
!	
	g = 1.0d0
	do while (g > 1d-10 .or. g < -1d-10)  
!
	C=0.5*(a+b) !! why is the do while block not indented? Also, use 0.5_dp instead of 0.5
	etastart = 0.0d0 !! this group of statements is repeated below. It could be made 
                     !! an internal subroutine
	etastop = 4.8d0 ! eta at infinity
	deta = 0.01d0
	y(1) = 0.0d0 ! Initial Condition f(0) = 0
	y(2) = 0.0d0 ! Initial Condition f'(0) = 0
	y(3) = c ! Initial Condition guess for f''(0)
	eta = etastart
!
		do while(eta < etastop)
			call Runge ()
			if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
			end if !! should there be a statement before this?
		end do
!
	g = y(2)-1d0
!
		if (g > 0d0)then !! put a space before then for readability
			a = c
			else !! should have same indentation as if
			b = c
		end if
	end do
	write (6,*)''
	Write (6,*) 'The value for f"(0) that satisfies fprime(inf) = 1'
	Write (6,'(a,f10.5)') 'f"(0) =' ,c
	Write (6,*)''
!
! Now calculate the solutions to the ODE's
!
	etastart = 0.0d0
	etastop = 4.8d0 ! eta at infinity
	deta = 0.01d0
	y(1) = 0.0d0 ! Initial Condition f(0) = 0
	y(2) = 0.0d0 ! Initial Condition f'(0) = 0
	y(3) = c ! Initial Condition f"(0) = c
	eta = etastart
	Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
	do while(eta < etastop)
			call Runge () !! should have same indentation as line below
		if (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
			Write (6, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
		end if
	end do
!
!	Program subroutines and sub functions 
!
Contains
	Subroutine Runge ()
		k1 = deta*f(eta,y) !Runge-cutta 4th order formulas
		k2 = deta*f(eta+deta/2, y+k1/2)
		k3 = deta*f(eta+deta/2, y+k2/2)
		k4 = deta*f(eta+deta, y+k3)
		y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
		eta = eta+deta
	end subroutine
!
	Function f(eta,y)
		real(dp), intent(in) :: eta, y(n)
		real (dp) :: f(n)
		f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
	end Function
!
end Program AE6009

```

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [November 25, 2021, 1:33pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/27 "2021-11-25T13:33:40Z")

</div>

Beliavsky,  
Thank you very much for your comments. I’ll work on them all.

Is there any way I can reduce or consolidate this output into perhaps one line?

```auto
    Write (output_unit,*)''
	Write (output_unit,*) "The value for f''(0) that satisfies f'(inf) = 1 is:"
	Write (output_unit,'(a,f10.5)') ' f"(0) =' ,c
	Write (output_unit,*)''

```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [November 25, 2021, 1:47pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/28 "2021-11-25T13:47:27Z")

</div>

> [@RickP330](#):
>
> Is there any way I can reduce or consolidate this output into perhaps one line?
> 
> ```auto
> Write (output_unit,*)''
> Write (output_unit,*) "The value for f''(0) that satisfies f'(inf) = 1 is:"
> Write (output_unit,'(a,f10.5)') ' f"(0) =' ,c
> Write (output_unit,*)''
> 
> ```

Yes, use the slash `/` edit descriptor, described for example [here](https://www.obliquity.com/computer/fortran/format.html). For example, first two lines above can be written

` Write (output_unit,"(/,a)") "The value for f''(0) that satisfies f'(inf) = 1 is:"`

Fortran formats have lots of obscure but useful features. I wonder what the best guide to learning them is.

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [November 25, 2021, 4:06pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/29 "2021-11-25T16:06:29Z")

</div>

Hi Guys,  
Okay, some typo’s aside, I think it’s done. Any other recommendations to tighten it up? I did notice that if I enter guesses too high I get NaN errors in the run. I am not sure how to avoid that.  
RP

```auto
Program AE6009
	use iso_fortran_env !To specify common output_unit
	implicit none !To declare all vars explicitly
!
	integer, parameter :: dp = kind(1d0)
	integer, parameter :: n = 3 !Order of ODE's
	real(dp) :: a, b, c, g !Variables for shooting method
	real(dp) :: eta, deta, etastart, etastop !Variables for Runge-cutta
	real(dp), dimension(n) :: y, k1, k2, k3, k4
	integer :: i
!	
	Write (output_unit,*)''
	Write (output_unit,*) "Enter a guess for f''(0) > 0"
	Read (*,*) a
	Write (output_unit,*) "Enter a guess for f''(0) < 0"
	Read (*,*) b
!	
!	Test for valid input variables
!
	If (a < 0 .or. b<0) then
	Write (output_unit,*) 'Negative values for guesses not permitted'
	Stop
	End if
!
	CAll Initialize_Variables()
	y(3) = a
		Do while(eta < etastop)
		Call Runge ()
		End do
			If (Y(2) < 0d0) then
			Write (output_unit,*) "Bad guess for f''(0) > 0"
			Write (output_unit,'(a,f10.5)') ' f"(0) = ' ,y(2)
			Stop
			End if
!
	Call Initialize_Variables()
	Y(3) = b
		Do while(eta < etastop)
		Call Runge ()
		End do
			If (Y(2) > 0d0) then
			Write (output_unit,*) "Bad guess for f''(0) < 0"
			Write (output_unit,'(a,f10.5)') ' f"(0) = ' ,y(2)
			Stop
			End if
!
	g = 1.0d0 !We need some IC for g
	i = 0
		Do while (g > 1d-10 .or. g < -1d-10)
		i = i+1
		C = 0.5d0 * (a + b)
		Call Initialize_Variables()
			Do while(eta < etastop)
			Call Runge ()
				If (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-10) then
				End if
			End do
		g = y(2)-1d0
			If (g > 0d0)then
			a = c
			Else
			b = c
			End if
				If (i > 1d4) then
				Write (output_unit,*) "Input guesses too far apart, infinite loop aborted"
				Stop
				End if
		End do
	Write (output_unit,"(/a)") "The value for f''(0) that satisfies f'(inf) = 1 is:"
	Write (output_unit,'(a,f10.5/a)') ' f"(0) =' ,c
!
! Now calculate the solutions to the ODE's
!
	Call Initialize_Variables()
	Write (output_unit, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
!
		Do while(eta < etastop)
		Call Runge ()
			If (abs(nint(10d0*eta)-(10d0*eta)) <= 1d-12) then
			Write (output_unit, '(a,f4.1,a,3f10.5)') 'y(',eta,') = ',y
			End if
		End do
!
!	Program subroutines and sub functions 
!
Contains
	Subroutine Initialize_Variables()
		etastart = 0.0d0
		etastop = 4.8d0 ! eta at infinity
		deta = 0.01d0
		y(1) = 0.0d0 ! Initial Condition f(0) = 0
		y(2) = 0.0d0 ! Initial Condition f'(0) = 0
		y(3) = c ! Initial Condition f"(0) = c
		eta = etastart
	End subroutine
!
	Subroutine Runge ()
		k1 = deta*f(eta,y) !Runge-cutta 4th order formulas
		k2 = deta*f(eta+deta/2, y+k1/2)
		k3 = deta*f(eta+deta/2, y+k2/2)
		k4 = deta*f(eta+deta, y+k3)
		y = y+((1d0/6)*(k1+(2*k2)+(2*k3)+k4))
		eta = eta+deta
	End subroutine
!
	Function f(eta,y)
		real(dp), intent(in) :: eta, y(n)
		real (dp) :: f(n)
		f = [y(2), y(3), y(2)*y(2)-1d0-y(1)*y(3)] ! Hiemenz Equation
	End Function
!
End Program AE6009

```

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [November 25, 2021, 6:06pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/30 "2021-11-25T18:06:18Z")

</div>

Id like to output the data to a CSV file or something I can import into a charting software. Is there a simple way to do this?  
Rick

---

<div class="post-metadata">

**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 26, 2021, 3:07am UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/31 "2021-11-26T03:07:56Z")

</div>

Something similar to

```auto
write (unit, 100) eta, y
100 format(F4.1,3(',',F10.5))

```

should suffice, assuming that you take care of opening the file before writing to it, and closing the file after all the writing has been done.

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [December 1, 2021, 12:29pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/32 "2021-12-01T12:29:59Z")

</div>

Hello All,  
Just to update you, I finished the program and submitted it. It works great! With all of your help it really made it much easier for me. Thank you.  
Only thing I didn’t resolve which was not necessary, if entered a value for f"(0) which was too large I got NaN errors in the output file. I assume those are negative infinity errors due to DP numbers getting so large, but not really sure. I just instruct the user to input closer values.  
If you have any comments on the NaN errors I would appreciate it.  
If there are any questions about the program please let me know. FORTRAN is a very powerful tool, I’m glad I was forced to revisit it. Looking forward to another project I can use it on…  
Kind regards,  
Rick

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [December 1, 2021, 12:43pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/33 "2021-12-01T12:43:49Z")

</div>

> [@RickP330](#):
>
> FORTRAN is a very powerful tool, I’m glad I was forced to revisit it. Looking forward to another project I can use it on…

Glad it worked for you. Note that since the 1990 standard, the official spelling is Fortran, with only the first letter capitalized. Wikipedia [says](https://en.wikipedia.org/wiki/Fortran)

> The names of earlier versions of the language through FORTRAN 77 were conventionally spelled in all-uppercase (FORTRAN 77 was the last version in which the Fortran character set included only uppercase letters. The official language standards for Fortran have referred to the language as “Fortran” with (rather than “FORTRAN” in all-uppercase) since Fortran 90.

---

<div class="post-metadata">

**Author:** ![RickP330](https://avatars.discourse-cdn.com/v4/letter/r/ecccb3/32.png) [@RickP330](https://fortran-lang.discourse.group/u/RickP330)\
**Post date:** [December 1, 2021, 6:30pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/34 "2021-12-01T18:30:31Z")

</div>

Hi Beliavsky,  
Thanks, as if I didn’t feel old enough already. 😵 The “kids” in my graduate class weren’t even born yet when I graduated with my UG. It’s funny to hear them all complain about having to program in a machine type language like this. None of them get it. I couldn’t even imagine them having to run a program from a terminal prompt instead of icon on the desk top.  
It was FORTRAN 77 I studied in 1989.  
Thanks for the help.  
Rick

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [December 1, 2021, 7:06pm UTC](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275/35 "2021-12-01T19:06:50Z")

</div>

> [@RickP330](#):
>
> … I couldn’t even imagine them having to run a program from a terminal prompt instead of icon on the desk top …

@RickP330 , are the other students in your class heavily into non-STEM pursuits for some reason? Or at commercial companies influenced by Windows OS or something?

Because otherwise with the increasing popularity of interpreted and interactive computing environments (Python, R. etc.) with the so-called [**REPL**](https://en.wikipedia.org/wiki/Read%E2%80%93eval%E2%80%93print_loop) mode and also Linux (so much work from command line), what people are doing now is conceptually not all that different from what I recently saw in a video at the enterprise where I work of late 1970s of employees doing calculations using various job scripts on terminals connected to an enterprise mainframe system.

[Previous page](https://fortran-lang.discourse.group/t/program-to-numerically-solve-3-odes/2275.md?page=1)