Hello everyone!
Please, somebody know how could I find the coincident point between the curves below?
I thought using two tables and use “if” in case the point is equal, but I don’t know how to do this and if this is the best way.
I need to find the value in X for the point in red.
Thank you!
If f(x) is the green function and g(x) the red one, what you want is solving f(x)-g(x)=0.
See Root-finding algorithms - Wikipedia
Note that equality will probably not be obtained perfectly for numerical reasons, but you can study \lvert f(x)-g(x) \rvert < \epsilon
Thank you, @vmagnin, for the answer. I tried to use root method, but the function is to complex and I didn’t reach a result for that. Because this I’m trying to do in other way.
E = 210000
n = 0.14
A = input by user
x=y/E+(y/K)^1/n
x=A^2/(E*y)
I tried to use www.wolframalpha.com to help me, but it isn’t work 
because this I was thinkin in create two array and find a way to compare them and find the point.
program intersection
implicit none
real, allocatable :: x_points(:)
real, allocatable :: f_vals(:)
real, allocatable :: g_vals(:)
integer :: intersection_index
real :: intersection_x, intersection_y
! ...
x_points = ...
f_vals = f(x_points)
g_vals = g(x_points)
intersection_index = minloc(abs(f_vals - g_vals))
intersection_x = x_points(intersection_index)
intersection_y = f_vals(intersection_index)
print *, intersection_x, intersection_y
end program
Where f and g are elemental functions of the form you show. That should get you pretty close. If you want to get a bit more fancy you could do some interpolation, sample more points in that range, etc. but hopefully you get the idea.
Your problem reduces to finding the root of a nonlinear equation in one variable. While it may be conceptually helpful to think of this problem as finding the intersection of two curves, you need not hamstring the solution process by replacing the function by an array.
You did not specify the values of A and K. I assumed (arbitrarily) that A=100, K = 0.5, and solved the problem using an old version of Maple:
E:=2.1e5: n:=0.14: K:=0.5: A:=100: fsolve(y/E+(y/K)^(1.0/n)-A*A/(E*y),y=0.4);
to obtain the result y = 0.3745933332.
Note that I used an initial guess of 0.4 for the solution. Some trial and error may be needed to find an initial guess that is good enough for a solution to be obtained. Depending on the values of E, n, K and A, there may be multiple solutions for y, or there may be no real solutions.
Thank you for all answers! @everythingfunctional I’m trying this way, but I’m having error in minloc syntax, could you help me, please? Following the code:
program Neuber_test2
implicit none
integer :: i
real :: x (25) = (/(i, i=1,25, 1)/)
real :: neuber (25), ramberg(25), resultado(25), raiz
do i=1, 25
neuber(i)=x(i)**2
end do
do i=1, 25
ramberg(i)=-x(i)**2+200
end do
do i=1, 25
resultado(i)=abs(neuber(i)-ramberg(i))
end do
raiz=minloc(resultado(i))
print*, raiz
end program Neuber_test2
Error: ‘array’ argument of ‘minloc’ intrinsic must be an array|
Thank you!
The error message is explicit: minloc is looking for the minimum element of an array. resultado(i) is scalar, not an array. You should write raiz=minloc(resultado) (or actually raiz=minloc(resultado,dim=1) if you want minloc to directly return a scalar)
Yep. I forgot about that. Always be wary of forum code. The compiler in my wetware isn’t perfect. 
I will note, Fortran is an array oriented language. All those loops are unnecessary (and try not to use magic numbers). I.e.
program Neuber_test2
implicit none
integer :: i
integer, parameter :: n == 25
real, parameter :: x (*) = [(i, i=1,n)]
real :: neuber(n), ramberg(n), resultado(n), raiz
neuber=x**2
ramberg=-x**2+200
resultado=abs(neuber-ramberg)
raiz=minloc(resultado, dim=1)
print*, raiz
end program Neuber_test2
While it’s not your actual problem, I did a subroutine to find crossings between two lines that don’t share the same x values. Maybe it can be helpful for someone looking at this topic for that case and can find it helpful
type :: cross
real(pr) :: x
real(pr) :: y
integer :: i
integer :: j
end type cross
subroutine find_cross(x_values1, x_values2, y_values1, y_values2, crossings)
!! Find crossings between two given lines
!!
!! Returns an array of crossigns, containings the crosses found. Each row
!! contains the data from each found cross
!!
!! | --------| ------- | ---------------- | ----------------- |
!! | x_cross | y_cross | first_line_index | second_line_index |
!! | --------| ------- | ---------------- | ----------------- |
!!
real(pr), allocatable, intent(in) :: x_values1(:) !! First line x values
real(pr), allocatable, intent(in) :: x_values2(:) !! Second line x values
real(pr), allocatable, intent(in) :: y_values1(:) !! First line y values
real(pr), allocatable, intent(in) :: y_values2(:) !! Second line y values
type(cross), allocatable :: crossings(:) !! Array of crossings
type(cross) :: current_cross
real(pr) :: x11, x12, x21, x22, y11, y12, y21, y22
real(pr) :: x_cross, y_cross, m1, b1, m2, b2, xlow, xup, ylow, yup
real(pr), dimension(2) :: xpair_1, xpair_2, ypair_1, ypair_2
integer :: i, j, n
if (allocated(crossings)) then
deallocate (crossings)
end if
allocate (crossings(0))
n = 0
do i = 2, size(x_values1)
xpair_1 = x_values1(i - 1:i)
ypair_1 = y_values1(i - 1:i)
x11 = xpair_1(1)
x12 = xpair_1(2)
y11 = ypair_1(1)
y12 = ypair_1(2)
m1 = (y12 - y11)/(x12 - x11)
b1 = y11 - m1*x11
do j = 2, size(x_values2)
xpair_2 = x_values2(j - 1:j)
ypair_2 = y_values2(j - 1:j)
x21 = xpair_2(1)
x22 = xpair_2(2)
y21 = ypair_2(1)
y22 = ypair_2(2)
m2 = (y22 - y21)/(x22 - x21)
b2 = y21 - m2*x21
x_cross = (b1 - b2)/(m2 - m1)
y_cross = m1*x_cross + b1
xlow = max(minval(xpair_1), minval(xpair_2))
xup = min(maxval(xpair_1), maxval(xpair_2))
ylow = max(minval(ypair_1), minval(ypair_2))
yup = min(maxval(ypair_1), maxval(ypair_2))
if ( &
(xlow <= x_cross) .and. (x_cross <= xup) .and. &
(ylow <= y_cross) .and. (y_cross <= yup) &
) then
print *, "CROSS:", i, j, x_cross, y_cross
if ((abs(x_cross - crossings(n)%x) < 0.1) .and. &
(abs(y_cross - crossings(n)%y) < 0.1)) then
print *, "CROSS: Repeated cross, skipping..."
cycle
end if
current_cross = cross(x_cross, y_cross, i, j)
n = n + 1
crossings = [crossings, current_cross]
end if
end do
end do
end subroutine find_cross
An alternative where x(i) does not have to be i.
This could be used to revise x range, where x = 10 is not an exact solution
program Neuber_test2
implicit none
integer :: i
real :: x (25) = (/(i, i=1,25, 1)/)
real :: neuber (25), ramberg(25), resultado(25), raiz
do i = 1, 25
x(i) = 8.5 + i/10.
neuber(i) = x(i)**2
ramberg(i) = -x(i)**2+200
resultado(i) = neuber(i)-ramberg(i)
write (*,*) i, x(i), resultado(i)
end do
i = minloc ( abs(resultado) ) ! stupid this is non-standard
raiz = resultado(i)
print*, "min at",i, x(i), raiz
end program Neuber_test2
For a more complex function where minval /= 0, you could look at repeating for new values of X, around the previous closest value.
Perhaps the array size does not have to be 25.
I think this is a good example of where automatic differentiation can be really useful since you only have to write your functions - the derivatives are handled for you:
program intersection
! Only need first derivatives from the autodiff library
use AVD, T1=>avd_d1
implicit none
! Equation parameters
real(8), parameter :: A = 100.0d0
real(8), parameter :: K = 0.5d0
integer, parameter :: E = 210000
! Solver parameters
real(8), parameter :: eps = 1.0d-8
integer, parameter :: max_iterations = 20
! Local variables
real(8) :: x
type(T1) :: r
logical :: solved
integer :: counter
! -------------------------------------------------------------------------
! Initialize autodiff
call init
! Calculate solution using NR (using associate to enhance readability)
x = 1.0d0
solved = .false.
associate (fx => r%v, dfx => r%d1)
newton_raphson: do counter=1,max_iterations
r = f(x)
if (abs(fx) < eps) then
solved = .true.
exit newton_raphson
end if
x = x - fx/dfx
end do newton_raphson
end associate
if (solved) then
write(*,'(a,f10.8)') 'Result = ',x
else
write(*,'(a)') 'FAILED'
end if
stop
contains
! Wrapper to simplify the NR code
type(T1) function f(x)
real(8), intent(in) :: x
type(T1) :: z
z = T1(x,1)
f = f1(z) - f2(z)
end function f
! x=y/E+(y/K)^1/n
type(T1) function f1(x)
type(T1), intent(in) :: x
real(8), parameter :: n = 0.14d0
f1 = x/E + (x/K)**(1/n)
end function f1
! x=A^2/(E*y)
type(T1) function f2(x)
type(T1), intent(in) :: x
f2 = A**2/(E*x)
end function f2
end program intersection
Here I’ve used my own AD library but no doubt all the others offer a similar interface. So the NR loop is written as in the text books and your functions are as you write them.
I chose the same parameters as in a previous reply with the output:
Result = 0.37459333