How can I solve this equation which has condition?

Help Homework
21

Thank @Arjen so much :heart_eyes:

22

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

image
image810×546 10.3 KB

And here is my code, when i compile it seems to be right, as the picture below
image
image975×602 139 KB 

        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

image
image578×562 36.1 KB

I don’t know why. Please help me!

23

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
24

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
25

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
26

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?

27

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