Can anyone help me rectifying the error shown while executing this program?

I'm trying to run the program given as follows; I received the following error

Attempted to access Kvals(:,2); index out of bounds because numel(Kvals)=1.

Error in TestAll>Fie/mod_simp (line 174)|**

                kermat(:,col) = Kvals(:,col)*wt(col);

Error in TestAll>Fie (line 109)

                [told,solnold] = mod_simp(n);

Error in TestAll>TestSmooth (line 23)

                [soln,errest,cond] = 
                Fie(lambda,a,b,behavior,@kernel,@RHS,AbsTol,RelTol);

Error in TestAll (line 4)

                TestSmooth(1,0,1)

       %==================The program is as follows;

       function TestAll

       disp('SMOOTH KERNEL')

       TestSmooth(1,0,1)

       disp(' ')

       disp('********************************************************************')

       disp(' ')

       end

       %==========================================================================

       function TestSmooth(lambda,a,b,AbsTol,RelTol)

       behavior =1;

       if nargin < 5 isempty(AbsTol)

       AbsTol = 1e-6;

       end

       if nargin < 6 isempty(RelTol)

       RelTol = 1e-3;

       end

       % Compute and plot the solution:

       [soln,errest,cond] = Fie(lambda,a,b,behavior,@kernel,@RHS,AbsTol,RelTol);

       t = soln.s; sol = soln.x;

       nfinal = length(t) - 1;

       % Interpolate the solution for assessing the error and if necessary,

       % use to get a smooth graph.

       tint = linspace(a,b,150);

       xint = ntrpFie(soln,tint);

       if nfinal < 150

       plot(tint,xint)

       else

       plot(t,sol)

       end

       % Study the error and condition of the integral equation.

       true_error = norm(true_soln(t) - sol,inf);

       max_error = norm(true_soln(tint)-xint,inf);

       disp(' ')

       disp(['Condition number = ',num2str(cond)])

       disp(' ')

       fprintf('Approximate bound on error at nodes = %6.1e\n',errest)

       fprintf('Actual error at nodes = %6.1e\n',true_error)

       fprintf('\nActual error at 150 equally spaced points = %6.1e\n',max_error)

       function kst = kernel(s,t)

       nn=16;

       height=20e-6; % half-height of the plate

       aa=10e-6;

       F3=@(ss)(1-exp(-2.*(ss./aa).*height))./((1+exp(-2.*
       (ss./aa).*height)));%F3(s/aa)

       sdata=0:0.1:1;tdata=0:0.1:1;

       for ii=1:length(sdata);

       s=sdata(ii);

       for jj=1:length(tdata);

       t=tdata(jj);

       f21=@(ss)ss.*(F3(ss)-1).*besselj(0,ss.*s.^2).*besselj(0,ss.*t.^2);

       f2=@(ss)-2.*t.*s.*t.*f21(ss);

       I(ii,jj)=galag(f2,nn);

       end

       end

       J=sum(I); % kernel

       kst=sum(J);

       function ts = true_soln(t)

       ts = 0;

       end

       function rs = RHS(s)

       rs= s ;

       end

       end

       function s = galag(func,n)

       if (n==2)|(n==4)|(n==8)|n(n==16)

       c=zeros(16,3); t=zeros(16,3);

       c(1:2,1)=[1.533326033   4.45095733]; % Cs are the weights

       c(1:4,2)=[0.832739124   2.048102438 3.631146306 6.487145084];

       c(1:8,3)=[0.43772341    1.033869348 1.669709766 2.376924702 3.208540913 
       4.268575511 5.818083369 8.906226215];

       c(:,4)=[0.225036315 0.525836053 0.831961392 1.146099241 1.471751317 
       1.813134687 2.17551752 2.56576275   2.993215086 3.471234483 4.020044086 
       4.672516608 5.487420658 6.585361233 8.276357984 11.82427755];

       t(1:2,1)=[0.585786438   3.414213562]; % ts are the nodes

       t(1:4,2)=[0.32254769    1.745761101 4.536620297 9.395070912];

       t(1:8,3)=[0.170279632   0.903701777 2.25108663  4.26670017  7.045905402 
       10.75851601 15.74067864 22.86313174];

       t(:,4)=[0.08764941  0.462696329 1.141057775 2.129283645 3.437086634 
       5.078018615 7.070338535 9.438314336 12.21422337 15.44152737 19.18015686 
       23.51590569 28.57872974 34.5833987  41.94045265 51.70116034];

       j=1;

        while j<=3

        if 2^j==n; break

        else

        j=j+1;

        end

        end

        s=0;

        for k=1:n

        x=t(k,j);y=feval(func,x);

        s=s+c(k,j)*y;

        end

        end

        end

        end

        %==========================================================================

        function [sol,errest,cond] = 
        Fie(lambda,a,b,behavior,kernel,RHS,AbsTol,RelTol)

        % Create and initialize solution structure sol.

        n = 2^3;

        [told,solnold] = mod_simp(n);

        if behavior == 0

        sol.s = (1 - told) ./ told;

        else

        sol.s = told;

        end sol.x = solnold;

        sol.lambda = lambda;

        sol.kernel = kernel;

        sol.behavior = behavior;

        if behavior == 4

        sol.alpha = alpha;

        end

        sol.RHS = RHS;

        max_m = 9;

        for m = 4:max_m

        n = 2^m;
        [t,soln] = mod_simp(n);   % Solve with n subdivisions of [a,b].

        % Calculate the norm of the difference between the current numerical
        % solution and the preceding one. In several cases we must interpolate
        % to get approximations on the same mesh.

        delta = norm(soln(1:2:n+1) - solnold,inf);

        % Update solution structure.

        sol.s = t;

        sol.x = soln;

        if m > 4

        % Estimate contraction rate and use to predict error.  
        % Be conservative and do not use if not contracting.

        rate = min(0.5,max(delta/deltaold,0.0625));

        end

        errest = (rate/(1 - rate))*delta; 

        if errest < (atol + rtol*norm(soln,inf))

        break;

        end  

        deltaold = delta;

        if m == max_m

        warning('Fie:Failure','Failed to pass the error test.')

        end

        end

        % Inexpensive approximation of condition of matrix in 1-norm.

        if nargout == 3

        cond = condest(kermat);

        end

        function [t,soln] = mod_simp(n)

        % Solve with n subdivisions of [a,b] using a quadrature rule. Simpson's

        % rule is used for smooth kernels.

        h = (b - a)/n;

        t = linspace(a,b,n+1)';

        % Weights for Simpson's rule.  Assumes n is even.

        wt = ones(1,n+1); wt(3:2:n-1) = 2; wt(2:2:n) = 4;

        wt = (h/3)*wt;

        [S,T] = meshgrid(t,t);

        rhs = RHS;

        Kvals = kernel(S,T)';

        Msize = n + 1;

        kermat = zeros(Msize);

        for col = 1:n+1

            kermat(:,col) = Kvals(:,col)*wt(col);

        end

        for j = 1:Msize

        kermat(j,j) = kermat(j,j) - lambda;

        end

        soln = kermat \ (-rhs);

        end % mod_simp

        end

Thanks in advance

Answer 1

[Jan]

The error message is clear: After

Kvals = kernel(S,T)';

Kvals contains one column only. Why do you assume that it has n+1 columns?

Did you use the debugger already? Type this in the command window:

dbstop if error

Now Matlab stops, when the error occurres and you can check the values of the locally used variables.

[Richardson Joseph]

Thanks for your answer. Kvals works well for a simpler case e.g. kst = (s.^3).*t.^4. However, for real case, the kernel is an infinite integral which I tried to solve using guass laguerre method.

[Richardson Joseph]

Now I have identified the error and the code is working well. Thanks Jan for the help.