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

integral equationsimpson's method

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

Best Answer

  • 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.