Help needed to calculate Jacobian by forward difference method

Help
1

I am trying to calculate the Jacobian of the two non linear equations using forward diff method. Here in this code I written a subroutine for the Jacobian. and function for the linear equations.
and main program to write result. when I execute the program I am getting wrong values for the function and Jacobian. I request to explain the code why its happening.
will attach the code below

real function Myfun(x, n, f)

implicit none
integer :: n

real :: x(n)
real :: f(n)
f(1) = x(1)**2 + x(2) - 11
f(2) = x(1) + x(2)**2 - 7

end function Myfun

program mainMNewton

implicit none
real, external :: Myfun
integer :: n
parameter(n=2)
integer :: np
parameter (np=4)
real, dimension(np) :: df
real, dimension(n) :: x


x = reshape((/3.2, -1.7/), shape(x))

call FdJac(n, x, Myfun, np, df)

write(*,*) "Jac:", df

end program

subroutine FdJac(n, x, fvec, np, df)
    integer :: n, np ,Nmax
    real :: df(np,np),fvec(n), x(n), EPS
    parameter (Nmax = 40, EPS = 1.e-4)


    integer :: i, j
    real :: h, temp, f(Nmax)



    do j = 1, n
        temp = x(j)
        h = EPS * abs(temp)
        if (h .eq. 0) h = EPS
        x(j) = temp + h
        h = x(j) - temp
        call Myfun(n, x, f)  ! call for the function vector
        x(j) = temp
        do i = 1, n
            df(i,j) = (f(i)-fvec(i)) / h
            write(*,*) "f= ", f(i), "fvec", fvec(i)
        end do
    end do
    return
end subroutine
2

I don’t know what’s the range of your variable, but this line seems wrong: I think you want to have an absolute threshold only if the relative value falls BELOW what can be represented on your machine:

if (h==0) h = EPS

In other words,

h = max(EPS, EPS*abs(temp))

Another thing I’ve noticed is that you never initialize the unperturbed solution, fvec: you should place an initial call at the beginning of FdJac:


! Initialize unperturbed vector
call Myfun(n,x,fvec)

do j=1,n
...
end do
3

its if(h .eq. 0) h = EPS. In call Myfun how can i include do loop if I have multivariable equation. can you brief me.

4

when include call Myfun(n,x,fvec) at the bigining of FdJac before
do j ) 1,n
…

the programm is not even compiling

5

Why not use a library?:

1 Like