Does my method for applying the trapezoid rule to a double integral look correct?

Hello,

I am working on some MATLAB code which involves using the Trapezoid rule to solve a double integral, and then plot the result in 3D coordinates. I have gotten my algorithm to work for a simpler equation, but for this more complex one that I am working with, I am running into some issues.

The equation I’m trying to solve is as follows:

The solution, which I can solve directly to check my results, is:

My algorithm uses a series of four nested loops:

for i=1:xpoints%x iterations
   for j=1:tpoints%t iterations
       for k=1:tpoints%theta iterations
           for l=1:tpoints%tau iterations

You can see my full code at the end of the post. The innermost loop (l) is for solving the inner integral and iterates over tau. The next loop out is for solving the outer integral (theta from 0 to t) and iterates over theta. The next loop out iterates over t points, and the outermost loop over x points.

The trapezoid rule is stated this way:

(where h=x2-x1)

When I plot the exact solution to the integral, I get this:

But when I plot my trapezoid rule code, I get this:

So clearly, there is an issue in solving for u values. All values of taux and u come out as either zero or very small negative or positive numbers (less than 10^-3), so I expect a problem inside the innermost loop (l), but can’t seem to find one. Maybe there’s an issue in how I’ve organized the loops themselves?

Here is my code. If my explanation isn’t helpful in answering a question you have about it, please ask anything:

xpoints=21;%number of data points in x dimension
tpoints=51;%number of data points in t dimension
x = 0:(1/(xpoints-1)):1; %defines number of data points in x
t = 0:(1/(tpoints-1)):1; %defines numeber of data points in t
u=zeros(tpoints,xpoints); %defines size of dependant variable matrix
taux=zeros(tpoints,xpoints); %taux contains the values of the inner integral solved
h=1/(tpoints-1);%define step size
for i=1:xpoints%x iterations
   for j=1:tpoints%t iterations
       for k=1:tpoints%theta iterations
           for l=1:tpoints%tau iterations; here, t(j)=t, t(k)=theta, and 
t(l)=tau
               %evaluate values for the inner integral:
               if l==1%for the first point, use the right hand rule
                   taux(l,i)=1/2*h*(x(i)*(t(l)*t(k)-
t(l)^2)+x(i)*(t(l+1)*t(k)-t(l+1)^2));
               elseif l==tpoints%for the last point, use the left hand rule
                   taux(l,i)=taux(l-1,i)+1/2*h*(x(i)*(t(l-1)*t(k)-t(l-
1)^2)+x(i)*(t(l)*t(k)-t(l)^2));
               else%for all other points, use the trapezoid rule
                   taux(l,i)=taux(l-1,i)+1/2*h*((x(i)*(t(l-1)*t(k)-t(l-
1)^2)+x(i)*(t(l)*t(k)-t(l)^2))/2 ...
                           +(x(i)*(t(l)*t(k)-t(l)^2)+x(i)*(t(l+1)*t(k)-
t(l+1)^2))/2);
               end
           end
           %evaluate values for the second (outer) integral:
           if k==1
               u(k,i)=1/2*h*((t(j)-t(k))*taux(k)+(t(j)-t(k+1))*taux(k+1));
           elseif k==tpoints
               u(k,i)=u(k-1,i)+1/2*h*((t(j)-t(k-1))*taux(k-1)+(t(j)-
t(k))*taux(k));
           else
               u(k,i)=u(k-1,i)+1/2*h*(((t(j)-t(k-1))*taux(k-1)+(t(j)-
t(k))*taux(k))/2 ...
                    +((t(j)-t(k))*taux(k)+(t(j)-t(k+1))*taux(k+1))/2);
           end
       end
   end
end
surf(x,t,u)%plot results
title('Approximate Solution')
xlabel('x')
ylabel('t')
zlabel('u')

Does anybody have an idea of what the problem is?

[Roger Stafford]

You need to explain how the variable 'x' entered into your expression for 'u'. There was no 'x' anywhere in the original double integral, just tau's, theta's, and t's. Your expression for 'u' without the 'x' would be correct.