Attempted to access c(4) ,,, lsqnonlin

Question

Michelle 2014-3-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/123232-attempted-to-access-c-4-lsqnonlin

评论： Michelle 2014-3-26

采纳的回答： Alan Weiss

在 MATLAB Online 中打开

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