Attempted to access c(4) ,,, lsqnonlin

2 visualizaciones (últimos 30 días)
Michelle
Michelle el 26 de Mzo. de 2014
Comentada: Michelle el 26 de Mzo. de 2014
Hi i'm getting two errors
  • Attempted to access c(4); index out of bounds because numel(c)=3
  • Error in lsqnonlin Caused by: Failure in initial user-supplied objective function evaluation. LSQNONLIN cannot continue.
Any suggestions on how I can solve there errors ?
function V2mainprogram
clc;
clear;
format long;
fid=fopen('referencerawdata.txt','rt');
c=textscan(fid,'%f %f %f',1);
Ro=cell2mat(c(1,1));
Ri=cell2mat(c(1,2));
theta0=cell2mat(c(1,3));
theta0=theta0.*pi./180; %converting to radians
c=textscan(fid,'%f %f %f %f','Headerlines',1);
fclose(fid);
Pexp=cell2mat(c(1));
riexp=cell2mat(c(2));
lambdaexp=cell2mat(c(3));
Fexp=cell2mat(c(4));
options=optimset('largescale','off','MaxFunEvals',1e100,'TolFun',1e-30,'TolX',1e-30,'MaxIter',5e3,'Algorithm','levenberg-marquardt','Display','off');
w=[0.4679139345726910 0.3607615730481386 0.1713244923791704 0.4679139345726910 0.3607615730481386 0.1713244923791704];
ksi=[-0.2386191860831969 -0.6612093864662645 -0.9324695142031521 0.2386191860831969 0.6612093864662645 0.9324695142031521];
model='G';
iniconstants=[0.1 0.1 0.1];
RadLambexp=[riexp lambdaexp];
constants=lsqnonlin(@(x) errPF(x,RadLambexp,Pexp,Fexp,Ro,Ri,theta0,ksi,w,model),iniconstants,[],[],options);
fprintf('\n');
fprintf('G model constants\n');
fprintf('c1=%f\n',constants(1));
fprintf('c2=%f\n',constants(2));
fprintf('c3=%f\n',constants(3));
iniguess=[riexp lambdaexp];
RadLamb=lsqnonlin(@(x) errPF(constants,x,Pexp,Fexp,Ro,Ri,theta0,ksi,w,model),iniguess,[],[],options);
[Pressure,Force]=PF(constants,RadLamb,Ro,Ri,theta0,ksi,w,model); %uses optimized values to get theoretical pressure
RadLambG=RadLamb;
PressureG=Pressure;
ForceG=Force;
figure(1);
subplot(2,1,1);
plot(Pexp.*1000,riexp,('-r')); hold on;
plot(PressureFung.*1000,RadLambFung(:,1),'-k',PressureTH.*1000,RadLambTH(:,1),'--',PressureG.*1000,RadLambG(:,1),'--');
ylabel ('Inner radius [mm]');
xlabel('Inner pressure [KPa]');
subplot(2,1,2);
plot(Pexp.*1000,riexp,('-r')); hold on;
plot(PressureFung.*1000,RadLambFung(:,2),'-k',PressureTH.*1000,RadLambTH(:,2),'--',PressureG.*1000,RadLambG(:,2),'--');
ylabel ('Longitudinal stretch ratio');
xlabel('Inner pressure [KPa]');
end
%*****************************
function error=errPF(c,RadLamb,Pexp,Fexp,Ro,Ri,theta0,ksi,w,model)
[Pt,Ft]=PF(c,RadLamb,Ro,Ri,theta0,ksi,w,model);
errp=(Pt-Pexp).*Ri^2;
errf=(Ft-Fexp);
error=[errp;errf];
end
%***********************************
function[Pt,Ft]=PF(c,RadLamb,Ro,Ri,theta0,ksi,w,model)
ri=RadLamb(:,1);
lambda=RadLamb(:,2);
P=zeros(length(ri),1);
F=zeros(length(ri),1);
ro=sqrt(((Ro.^2-Ri.^2).*theta0)./(pi.*lambda)+ri.^2);
for i=1:length(ksi)
r=((ro+ri)./2+((ro-ri)./2).*ksi(i));%gd
R=sqrt(Ri.^2+(((pi.*lambda)./theta0).*(r.^2-ri.^2)));%gd
Ezz=0.5*((lambda.^2)-1);%gd
Ett=0.5*((((pi.*r)./(theta0.*R)).^2)-1);%gd
calcPartials=str2func(model);
[dWtt,dWzz]=calcPartials(c,Ett,Ezz);
P=P+w(i).*(((((pi.*r)./(theta0.*R)).^2).*dWtt)./r);%gd
F=F+w(i).*((2.*(lambda.^2).*dWzz)-(((pi.*r)./(theta0.*R)).^2).*dWtt).*r;%gd
end
Pt=((ro-ri)./2).*P;
Ft=pi.*((ro-ri)./2).*F;
end
%***************************************
function[dWtt,dWzz]=G(c,Ett,Ezz)
dWtt=1/2.*c(1).*((2.*c(2).*Ett)-(c(3).*(1./(((2.*Ett+1).^2).*(2.*Ezz+1)))).*(1./((2.*Ett+1).*(2.*Ezz+1)))).*(exp((c(2).*(Ett.^2))+(c(3).*(Ezz.^2))+((c(4)).*(1/4).*(((1./((2.*Ett+1).*(2.*Ezz+1))))-1.^2))));
dWzz=1/2.*c(1).*((2.*c(3).*Ezz)-(c(3).*(1./(((2.*Ezz+1).^2).*(2.*Ett+1)))).*(1./((2.*Ett+1).*(2.*Ezz+1)))).*(exp((c(2).*(Ett.^2))+(c(3).*(Ezz.^2))+((c(4)).*(1/4).*(((1./((2.*Ett+1).*(2.*Ezz+1))))-1.^2))));
end

Respuesta aceptada

Alan Weiss
Alan Weiss el 26 de Mzo. de 2014
Well, I see in
function[dWtt,dWzz]=G(c,Ett,Ezz)
lines with
...+((c(4)).*...
If c has only three components, then these lines will definitely thow an error.
Alan Weiss
MATLAB mathematical toolbox documentation

Más respuestas (0)

Categorías

Más información sobre Historical Contests en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by