Answers to exercise 1

1.
program meanvalue
! calculate the mean of three numbers
implicit none
real :: x,y,z,mean
write(*,*) "give three numbers"
read (*,*) x,y,z
mean = (x+y+z)/3.0
write(*,*) "mean is", mean
end program meanvalue	

2.
program quadratic
! solve the quadratic equation ax^2+bx+c=0
implicit none
real :: a,b,c, D, x1, x2, re, im
write(*,*) "give the coefficients a,b,c"
read (*,*) a,b,c

if (a==0) then      ! it is a linear equation
  if ((b==0) .and. (c==0)) then
     write (*,*) "any value satisfies the equation"
  else if ((b==0) .and. (c /= 0)) then
     write (*,*) "no solution"
  else
     x1 = -c/b
     write (*,*) "single root", x1
  end if
else                 ! a truly quadratic equation
  D=b**2-4*a*c       ! discriminant
  if (D>=0) then     ! roots are real
    x1 = (-b + sqrt(D))/(2*a)
    x2 = (-b - sqrt(D))/(2*a)
    write(*,*) "real roots", x1, x2
  else                ! complex roots
    re =  -b/(2*a)
    im = sqrt(-D)/(2*a)
    write (*,*) "complex roots", re, "+-i*", im
  end if
end if

end program quadratic



3.
program solve
! find the root of sin(x)-x+2=0
implicit none
real :: x0, x1
x0 = 1.5 ! guess initial value
x1 = asin(x0-2)
! iterate until the values do not change
do while (abs(x1-x0) > 1e-5)
x0 = x1
x1 = asin(x0-2)
write(6,*) x1
end do
write (6, *) x1, sin(x1)-x1+2
end program solve

output:

             NaN
             NaN             NaN

does not work!

try another version

program solve
! find the root of sin(x)-x+2=0
implicit none
real :: x0, x1
x0 = 1.5 ! guess initial value
x1 = sin(x0)+2
! iterate until the values do not change
do while (abs(x1-x0) > 1e-5)
x0 = x1
x1 = sin(x0)+2
write(6,*) x1
end do
write (6, *) x1, sin(x1)-x1+2
end program solve

output:

   2.1435995    
   2.8403850    
   2.2966738    
   2.7479172    
   ...

   2.5542006    
   2.5541921    
   2.5541921       7.15255737E-06

convergence is pretty slow



4.
program harmonic
! find a partial sum of the harmonic seris
implicit none
real :: term, sum=0.0, invsum=0.0
integer (selected_int_kind(8)) :: i, n
! just to be sure we can handle a large number of terms
write(6,*) "number of terms?"
read(5,*) n
do i=1,n
  sum = sum + 1.0/i ! note 1.0, not 1
end do
do i=n,1,-1
  invsum = invsum + 1.0/i 
end do
write (6, *) n, sum, invsum
end program harmonic


>./a.out
number of terms?
1
1 1.000000 1.000000

>./a.out
number of terms?
10
10 2.928968 2.928968

>./a.out
number of terms?
100
100 5.187378 5.1873770

> ./a.out
number of terms?
1000000
1000000 14.392652 15.40368

> ./a.out
number of terms?
10000000
10000000 15.40368 16.686031

> ./a.out
number of terms?
20000000
20000000 15.403683 17.390091    

The series is divergent. However, eventually terms
become so small that they will not affect the sum
that is several orders of magnitude bigger.