matlab program help looping over t and the indices of a matrix to fill in entries

so this program seems just to iterate over the first row i need it to iterate over the matrix, and since i couldnt get a true summation to work im adding matrices ad each iteration of t until 1500 so i can get a reasonable approximation here is the code

nx=21;                  %Number of steps in space(x)
ny=26;                  %Number of steps in space(y)       
niter=1500;             %Number of iterations 
dx=.1;                  %Width of space step(x)
dy=.1;                  %Width of space step(y)
x=0:dx:2;               %Range of x(0,2) and specifying the grid points
y=0:dy:2.5;             %Range of y(0,2.5) and specifying the grid points
t=1;                    %Initialize t
%Initial Conditions
p=zeros(ny,nx);         %Preallocating p
pn=zeros(ny,nx);        %Preallocating pn
ps=zeros(ny,nx);        %Preallocating ps
%Analytical solution matrix formation
%Explicit iterative scheme with indices and t 

for j=1:nx
    for k=1:ny
        for t=1:niter
            p(k,j)=(((4*1.35)/pi)*(sin(((2*t)-1)*pi*(j-1)/20) ...
               *sinh(((2*t)-1)*pi*(k-1)/20)/sinh(((2*t)-1)*pi*(20/25))));
            ps(k,j)=pn(k,j)+p(k,j);
            pn(k,j)=ps(k,j);
        end
    end
end
%Plotting the solution
contour(pn)

when i look at the variable pn or ps it has NaN values for everything except the first row

[per isakson]

Here are some links on debugging in Matlab

• Debugging in MATLAB
• Debugging session
• Debugging Approaches

Answer 1

[Roger Stafford]

When you get to the third row where k = 3 and where t has ascended to 1500, the value of

   sinh(((2*t)-1)*pi*(k-1)/20)

in the numerator and

   sinh((2*t-1)*pi*20/25)

in the denominator have both overflowed to plus infinity and their quotient would therefore be a NaN. To avoid this, you should use the identity

   sinh(A)/sinh(B) = (exp(A-B)-exp(-A-B))/(1-exp(-2*B))

so that applying this to

   sinh((2*t-1)*pi*(k-1)/20)/sinh((2*t-1)*pi*20/25)

will no longer produce a NaN. That will get you down to about the 18th row. After that even this quantity produces an overflow and you will thereafter get infinities in p, ps, and pn. No identity can prevent infinities in these last eight rows, but at least there will be no NaNs.

[Roger Stafford]

When I said "No identity can prevent infinities in these last eight rows, but at least there will be no NaNs" I forgotten that the sine factor will be changing the signs of the infinity quantities and this will also produce a NaN from the operation infinity minus infinity.

Answer 2

[Image Analyst]

Standard debugging procedure is to break up that ungodly long equation into smaller individual terms. Then check each term to see if it's valid. Step through the code with the debugger examining terms to see when/where/why one or more of the terms is nan.

Answer 3

[christopher zangler-scaduto]

would it help if i broke the trig terms into exponentials

Answer 4

[christopher zangler-scaduto]

so unfortunately non of the solutions so far have solved my problem can some one show me code that will work i was thinking mesh plot with vectors for j, and K and just hard coding a finite sum can someone show me syntax for that or any other method that will work

[Roger Stafford]

As I have already indicated to you, Christopher, by the time you have reached the 18th row or thereabouts, the ratio of the two hyperbolic sines will have exceeded the maximum value allowed for double precision floating point numbers. It is a fundamental difficulty of your algorithm. There is no remedy for this other than to use higher precision numbers, which can be done by using symbolic numbers of the Symbolic Toolbox, or by using for example John D'Errico's high precision package in the file exchange.