Solving 2D Integrals Using Fortran 95 and the Trapezoidal Method

[CMJ96]

Homework Statement

I'd like to solve the following integral using the 2D trapezoidal method in fortran 95
$$I=\int^{1.40406704}_{-1.40406704} \int^{x+1.40406704}_{x-1.40406704} exp(x+y) dy dx$$

Homework Equations

$$I= \frac{h}{2} \left(f_0+ 2 \sum_{i=1}^{n-1} f_i +f_n \right)+O(h^2) $$

The Attempt at a Solution

So I have wrote my program below, attempting to use a nested do loop to calculate the integral over the two dimensions. However I am not sure how to use the result of my do loop to arrive at a final answer, I'm also not sure if my loop is entirely correct, here is a copy of my script

program trap2d
implicit none
real :: hy, hx, c, d, integralx, integraly, func, I_numeric,x,y,integral
integer :: nx,ny,ix,iy

c=-(1.140406704)
d=(1.140406704)
ny=1000
nx=1000
hx=(d-c)/nx
integralx=0.0

do ix=1,nx-1
x=c+ix*hx
hy=((x+1.140406704)-(x-1.140406704))/ny
integraly=0.0

do iy=1,ny-1
y=x-1.140406704+iy*hy
integraly=integraly+hy*func(x,y)
end do

integralx=integralx+hx*integraly
end do

end program

function func(x,y)
implicit none
real :: func,x,y
func=exp(x+y)
end function func

 

[scottdave]

From what I follow in your code, it appears to be doing the integration sums, but I don't think it's doing trapezoidal. It appears to be just doing (left point) sums. I don't see it utilizing your relevant equation either. Can you go into a little more detail about your relevant equation, and how it's used to approximate the integral?

 

[tnich]

Your code looks pretty good. I think you have forgotten to add in the function values at the endpoints ($f_0/2$ and $f_n/2$) in both in the inner and outer do loops. Aside from that, I think you need a print statement to show the value of integralx before your end program statement.

 

[CMJ96]

Can you go into a little more detail about your relevant equation

So the relevant equation breaks the integral up into trapeziums, the area of a trapezium is the following $$A=h\frac{f_i+f_{i+1}}{2} $$and sums them up to find an approximation, the more strips the more accurate.
The expansion of the integral goes along the following lines,
$$I= \frac{hf_0}{2}+hf_1+hf_2+...hf_{n-1} +\frac{h}{2}f_n$$
Just to clarify this is the composite 2D trapezoid method

I think you have forgotten to add in the function values at the endpoints (f0/2f0/2f_0/2 and fn/2fn/2f_n/2) in both in the inner and outer do loops

For a two variable function how do I define the end points? in my case would it be $$f(1.40406704, 2(1.40406704)) ,f(1.40406704, 0), f(-1.40406704,0), f(-1.40406704, 2(-1.40406704)) $$ ?

 

[scottdave]

Now I see it, thanks. So as @tnich pointed out, you would need the f₀/2 and f_n/2 to make it a trapezoidal method.
Then you need to print the variable integralx

 

[CMJ96]

Now I see it, thanks. So as @tnich pointed out, you would need the f₀/2 and f_n/2 to make it a trapezoidal method.
Then you need to print the variable integralx

How exactly do I define $$f_n$$ and $$f_0$$ in 2D?

 

[DrClaude]

How exactly do I define $$f_n$$ and $$f_0$$ in 2D?

You are doing nested 1D integrals. At each iteration in ix, you need to add the case for iy = 0 and iy = ny. Then, in addition to the loop over ix, you need code for ix = 0 and ix = nx.

 

[CMJ96]

Would adding in func(x, x+1.140406704) and func(x,x-1.140406704) both multiplied by a factor of hy/2 for the y loop, and would func(1.140406704,y) and func(-1.140406704,y) multiplied by hx/2 work for the x loop?

 

[DrClaude]

Would adding in func(x, x+1.140406704) and func(x,x-1.140406704) both multiplied by a factor of hy/2 for the y loop, and would func(1.140406704,y) and func(-1.140406704,y) multiplied by hx/2 work for the x loop?

In the second case, you still have to loop over y.

The simplest way to do it (from a programmer's point of view) is to extend the loops to 0 to nx (and 0 to ny), and use an if statement to set the multiplication factor to hx or hx/2 (and correspondingly hy or hy/2) depending in the value of ix (iy).

 

[CMJ96]

Thank you so much guys, you've helped a lot, my integral now works to within 0.01 of the analytical solution. I appreciate the help so much :)