How can I solve this equation which has condition?

Hm, you do not define the type of the function. Can you post the complete program, as it is not clear from your message from what point in the code the error originates.

1 Like

As a meta comment on the appearance of your codes shared, I would like to suggest to use indentations when it comes to if-clauses, do-iterations, etc.

Contrasting to other languages (e.g., Python) indentation (probably most often) is functional irrelevant to modern Fortran. However, this form to structure code (cf. Arjen’s improvement of bisect.f90) visually into parts sharing a task/within a logical step helps to read the code. It only is a matter of complexity of the code written and time since initial design of the program that this equally may be helpful to you, too.

If the environment you use does not support this by default, I suggest to pass your code to a code formatter. An equivalent to yapf3 and black in Python is e.g., Peter Seewald’s fprettify, which recognizes a selection of these small units and indents them accordingly (this includes nested loops). While probably not designed to be a code linter (like pylint), loops not properly terminated are easily identified in fprettify's output, too.

1 Like

And here is my code:

program solvingpr1
implicit real*8(a-h,o-z)
integer, parameter:: n=1000
double precision ::x1,x2,esp
double precision:: roots(n)
integer i,nroots
x1=0.0
x2=20.0
eps=1.0e-7
call solving(f,x1,x2,eps,n,roots,nroots)
if (nroots==0) then
write(*,*) 'has no root' 
stop
end if
write (*,*) 'roots of this equation are'
do i=1,nroots,1
write(*,*) i, roots(i), f(roots(i))
end do
contains
function f(x)
implicit real*8(a-h,o-z)
double precision f,x
f=(-50/(1.0+exp(x-5.84)/0.5))+(12.33*(3.0-x**2/34.1056)+1.41/(x**2)
else
f=(-50/(1.0+exp(x-5.84)/0.5))+144/x**2+1.41/(x**2)
endif
end function f(x)

end program solvingpr1


!----------------------------------------------subroutine
subroutine solving(f,x1,x2,esp,n,roots,nroots)
implicit none
integer, parameter:: n=1000
double precision f,x1,x2,eps,roots(n)
double precision a,b,c,dx,root
integer n,i,j,nroots
interger, parameter::sk =200
dx = (x2-x1)/n
nroots=0
do j =1,n
a=x1+(j-1)*dx
b=a+dx
if (f(a)*f(b))>0) cycle
do i = 1, sk
c=(b+a)/2.0
if (f(c).f(a)<=0) then
b=c
else
a=c
end if
if (abs(a-b)<=eps) exit
end do
root=(a+b)/2.0
if (abs(f(root)<1.0) then
nroots=nroots+1
Roots(nroots)= root
end if
end do
end subroutine solving

I think I could solve this problem by the way if x=x1 we have f(x1)=a
then i subtract f(x1)-4.290 then I compare the result with a value which I can accept (eps = 1e-07). and then infer to the root :).
I repeat this for many time. Is it possible?

Thanks for posting that code, but:

  • Pick up on the explicit advice by @nbehrnd - it really really helps to understand the code!
  • Do not use “implicit real*8(a-h,o-z)” - it is a sure method to let simple bugs into your code. Use “implicit none” instead. It forces you to declare all variables, but it prevents stupid mistakes like typing “o” instead of “0”, simply because you did not declare “o” to be a variable. (And do not use “o” as a variable name, by the way)

Looking at you code: “interger” is not a Fortran keyword. And the function f in solving is not defined as a function but as a variable. Hence the call is inconsistent.

1 Like

thank you very much :slight_smile: for your nice comments. It really helps me!

It is good to compile with picky compiler options that point out problems in code that is still legal. With gfortran I suggest -Wall -Wextra -Wconversion-extra For a snippet of your code

implicit none
double precision :: x2,eps
x2 = 20.0
print*,x2
end

gfortran -Wconversion-extra says

xxdouble.f90:3:3:

    3 | x2=20.0
      |   1
Warning: Conversion from 'REAL(4)' to 'REAL(8)' at (1) [-Wconversion-extra]
xxdouble.f90:4:4:

    4 | eps=1.0e-7
      |    1
Warning: Conversion from 'REAL(4)' to 'REAL(8)' at (1) [-Wconversion-extra]

The Quickstart Fortran Tutorial, which I suggest reading, gives the example (slightly abbreviated)

program float
use, intrinsic :: iso_fortran_env, only: dp=>real64
implicit none
real(dp) :: float64
float64 = 1.0_dp ! Explicit suffix for literal constants
end program float

and gives the advice, which I endorse,

Always use a kind suffix for floating-point literal constants.

P.S. Not directed to OP. Should not gfortran with -Wall -Wextra include -Wconversion-extra ? Does FPM warn by default about double precision variables set to single precision literal constants?

2 Likes

@haidangwy ,

An inquiry first: are you a student? and is this a homework problem?

If yes, I suggest studying the numerical computing advancements related to the works of Dekker, Brent et al. on finding zeroes of functions c.f. https://blogs.mathworks.com/cleve/2015/10/26/zeroin-part-2-brents-version/#d354b12f-235d-474e-be05-d4dba9fa6fcc

2 Likes

Thanks Mr.Beliavsky. I’m a student. This is my assignment. I try to write a code for myself with Fortran.
Thank your for your link. I’m going to read this.

I want to thank everybody for helping me :wink:

It’s really interesting. When I solved some problem then new problems come. Today I’ve got a new problem that I want to take my roots from this equation. Let’s say I have 3 roots: 1,2,3. I need to assign them in a,b,c for the further calculation in this program.

a= root(1)
b= root(2)
c= root(3)

How can I do it?

Well, assuming your solution routine returns multiple roots in an array “root” (as was the case in your first post), then that is exactly the syntax you need.

2 Likes

amazing. I did it. I’m happy about it :slight_smile:

Thank @Arjen so much :heart_eyes:

Next, I’ve got trouble with solving integration.
I like to solving this equation, as the picture below


And here is my code, when i compile it seems to be right, as the picture below

        program integral_wb
	implicit none
        double precision f,x1,x2
        integer i,m
        real*8:: TP
        external f
        x1=0.02248
        x2=4.487679
        m=2
        do i=1,16
        call integral1D(f,x1,x2,TP,m)
        write (*,*) , real(TP)
        m=m*2
        end do
        end program integral_wb
        !use simpson rule 1D
        subroutine integral1D(f,x1,x2,TP,m)
        implicit none
        double precision f,x1,x2,TP,s
        double precision h,x
        integer m,i
        if (mod(m,2).ne.m) then 
        m=m+1
        end if
        s=0.0
        h=(x2-x1)/dfloat(m)
        do i=2,m-2,2
        x=x1+dfloat(i)*h
        s=s+2.0*f(x)+4.0*f(x+h)
        end do
        TP=(s+2.0*f(x)+4.0*f(x+h))*h/3.0
        return
        end subroutine integral1D
        
! my function
        function f(x)
        implicit none
        double precision f,x
        real::termf10,termf11
        termf10=3585.39622/197.33*3.0e23
        termf11=2.0*3585.39622/(197.33**2.0)
        f=termf10/sqrt(termf11*((-50.0d0/(1+(exp(x-5.373)/0.5)))+
     3 (1.44*50.0/5.373)*(3.0d0-(x/5.373)**2.0)+
     4 (197.33/(8.0*3585.39622*x**2.0))-4.290))
        end function f

The problem is NaN appearing


I don’t know why. Please help me!

Have you plotted the function you are trying to integrate over the domain of integration? Or at least printed out its values on a fine grid and examined them?

On an unrelated note, it’s good to compile with picky options, as I previously mentioned. You use the dfloat() extension, for which

gfortran -Wall -Wextra -std=f2018 dfloat.f90

on

print*,dfloat(1)
end

says

dfloat.f90:1:13:

    1 | print*,dfloat(1)
      |             1
Warning: The intrinsic 'dfloat' at (1) is not included in the selected standard but a GNU Fortran extension and 'dfloat' will be treated as if declared EXTERNAL.  Use an appropriate '-std='* option or define '-fall-intrinsics' to allow this intrinsic. [-Wintrinsics-std]
dfloat.f90:1:13: Warning: The intrinsic 'dfloat' at (1) is not included in the selected standard but a GNU Fortran extension and 'dfloat' will be treated as if declared EXTERNAL.  Use an appropriate '-std='* option or define '-fall-intrinsics' to allow this intrinsic. [-Wintrinsics-std]

You could use the dble intrinsic function instead, or preferably real with a kind argument.

2 Likes

How many times is this loop executed?

do i=2,m-2,2

What happens in the next statement when x is not defined?

TP=(s+2.0*f(x)+4.0*f(x+h))*h/3.0

1 Like

Just as notes aside:

  • Perhaps a symbolic integration of T were beneficial to the subsequent computation. I’m not sure if Fortran is suitable for this, nor if mathematicians would pick one among programming languages/programs like Maple, SymPy, Wolfram|Alpha.

  • The definitions term10, term11, f and f1 require already existing comprehension of the overall computation ahead. It is easier for reading the code (and later share/maintenance) if names of variables, functions, etc. are more self-explicatory.

  • Lines 69 till and including 71 in the screen photo look like an (overly long) definition of f1. Based on the advancement of the line count, however, this looks broken.

    For me, it worked easier to define short individual terms (with mnemotic names) computed separately from each other. And then to use these bits and bolts for the computation of the eventually interesting thing. As an illustration, I refer to a pattern similar to

    do while (delta > 0.005_rp)
       call compute_Vn(x, Vn)
       call compute_Vc(x, Vc)
       call compute_Vl(x, Vl)
    
       V = Vn + Vc + Vl
    
       ! update for the stop condition
       delta = V - 4.290_rp
    
       ! updated report to the CLI
       write (*, '(F6.4, 2F8.3, 3F12.3)') x, Vn, Vc, Vl, V, delta
    
       ! care for the advancement of x
       x = x + increment
    end do
    
    

    This allows an easier identification and correction of errors caused by a typo (or passing a division by zero) while gradually assembling the computation than working at once with a very long equation difficult to understand right at first glance. It still is possible to mute the intermediate reporting line as a comment, and to add one line to report only x, delta, and V

1 Like

Thank for it. I had repeated 18 times.
Wow, I forgot to consider it. If x make f(x) not define how can pass this?

I’m going to try your way. Thank you so much :heart_eyes: