How to verify in Fortran whether an iterative formula of a non-linear system will converge

Question

How to verify in Fortran whether an iterative formula of a non-linear system will converge to the root near (x,y)?

It was easy for a programming language which support symbolic computations. But how to do that in Fortran? Like getting partial derivative of the component functions and check whether they bounded near the root. But I couldn't do that in fortran or haven't the idea how to do that. It will be a great help for me if anyone give me some idea for the following non-linear system now or if possible for a general case.

I want to use Fixed point iteration method for this case

Main system:

x^2+y=7
x-y^2=4

Iterative form (given):

X(n+1)=\sqrt(7-Y(n)),
Y(n+1)=-\sqrt(X(n)-4),
(x0,y0)=(4.4,1.0)

---

Theorem (which I follow)

The issue is, I need to check the boundedness of the partial derivatives of \sqrt(7-Y) and -\sqrt(X-4) on some region around (x0,y0)=(4.4,1.0). I can write the partial derivative function in fortran but how to evaluate so many values and check it is bounded around the (4.4,1.0).

Update

One possibly right solution would be to get arrays of values around (4.4,1.0) like (4.4-h,1.0-h)*(4.4+h,1.0+h) and evaluate the defined partial derivative function and approximate their boundedness. I haven't encounter such problem in Fortran, so any suggestion on that also can help me a lot.

Answer 1

If you just want to check the boundedness of a function on a grid, you can do something like

program verify_convergence
  implicit none
  
  integer, parameter :: dp = selected_real_kind(15, 307)
  
  real(dp) :: centre_point(2)
  real(dp) :: step_size(2)
  integer :: no_steps(2)
  real(dp) :: point(2)
  real(dp) :: derivative(2)
  real(dp) :: threshold
  
  integer :: i,j
  real(dp) :: x,y
  
  ! Set fixed parameters
  centre_point = [4.4_dp, 1.0_dp]
  step_size = [0.1_dp, 0.1_dp]
  no_steps = [10, 10]
  threshold = 0.1_dp
  
  ! Loop over a 2-D grid of points around the centre point
  do i=-no_steps(1),no_steps(1)
    do j=-no_steps(2),no_steps(2)
      ! Generate the point, calculate the derivative at that point,
      !    and stop with a message if the derivative is not bounded.
      point = centre_point + [i*step_size(1), j*step_size(2)]
      derivative = calculate_derivative(point)
      if (any(abs(derivative)>threshold)) then
        write(*,*) 'Derivative not bounded'
        stop
      endif
    enddo
  enddo
  
  write(*,*) 'Derivative bounded'
contains
  
  ! Takes a co-ordinate, and returns the derivative.
  ! Replace this with whatever function you actually need.
  function calculate_derivative(point) result(output)
    real(dp), intent(in) :: point(2)
    real(dp) :: output(2)
    
    output = [sqrt(7-point(2)), -sqrt(point(1)-4)]
  end function
end program

I know the function calculate_derivative doesn't do what you want it to, but I'm not sure what function you actually want from your question. Just replace this function as required.

Answer 2

The main question is different: How can you calculate the solution of the mathematical problem without the help of any software? If you know that, we can program it in fortran or any language.

In particular, and assuming that n=0,1,2,3... to solve your problem you need to know X(0) and Y(0). With this you calculate

X(1)=\sqrt(7-Y(0)), Y(1)=-\sqrt(X(0)-4)

Now you know X(1) and Y(1), then you can calculate

X(2)=\sqrt(7-Y(1)), Y(2)=-\sqrt(X(1)-4)

etc.

If your system of equations converge to something, until some iterations (for example, n=10, 20, 100) you going the check that. But because the nature of fortran, it will not give you the solution in a symbolic way, it is not its goal.