Mostrar comentarios más antiguos
Hi, I have some problem using lsqcurvefit. It should work, as this is a simple case. I cannot find the mistake, so it hade been nice if somebody can see what i am doing wrong.
xdata is the time, for every 10 seconds. Nineos is the measured data, and the dimensions of xdata and Nineos is the same (3163x1)
xdata=[0:10:31620]';
Nineos=data((359:3521),6)*273.15; %(ydata)
param0=300;
options = optimset('Display','off','MaxIter',1e3,'MaxFunEvals',1e10,'TolFun',1e-6,'TolX',1e-6,'DiffMinChange',10);
[x,resnorm,residual,exitflag,output]=lsqcurvefit(@curvfit_U,param0,xdata,Nineos,0,700,options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function F=curvfit_U(param0,xdata)
U=param0(1);
load 2oct10
par = emulsion_parameters();
format('shortG')
x10=[350];
Nm_in=139204.4545;
Tr=data((359:3521),7)+273.15;
Tjin=data((359:3521),5)+273.15;
X=data((359:3521),2)/100;
Ts=10;
Nsample = length(Tr);
%%%%% ODE solver
x_ode = [];
t_ode = [];
for i=1:Nsample
if(i*Ts < 4840)
Emu=0; Ve2=0; m_em2_tot=0; Win=0; cin=0; Tin=0;
elseif (i*Ts > 4840)
Emu=492.57; Ve2=540; m_em2_tot=540000; Win=0.025796; cin=4.187; Tin=300;
end
t1 = Ts*(i-1);
t2 = Ts*i;
ts= 1;
tspan = t1 : ts : t2;
%X_in=X(i);
Tj_in=Tjin(i);
Tr_in =Tr(i);
[x1 time]= ode1(@estimation_U,tspan,x10,par,Emu,Ve2,m_em2_tot,Win,cin,Tin,Tr_in,Tj_in,U);
x10=x1(end,:);
x_ode = [x_ode ; x1];
t_ode = [t_ode time];
end
%t=t_ode';
toc
%time=[t(1:11:end)];
%state=x_ode;
state1=x_ode(1:11:end);
F=state1;
end
The last file is the model with the ODE equation:
function F =estimation_U(t,x,par,Emu,Ve2,m_em2_tot,Win,cin,Tin,Tr_in,Tj_in,U)
Tjin=Tj_in;
Tr=Tr_in;
%% Polynom fit to some data
Wc=((2.273*10^-13)*t^4)-((1.153*10^-8)*t^3)+((0.0001704)*t^2)-(0.6442*t)+6.244*10^4;
%% Part II. Extract present value of states Tjout=x(1); %% Part III. Intermediate calculations mj=1400000;
cpw=4.02*exp((1.99*10^-4)*Tr);
Tave=(Tjout+Tjin)/2;
%% Part III. Evaluate ODE set dx1=(((Wc*cpw*(Tjin-Tjout))+(par.Ua*par.Aa*(par.Ta-Tjout))+(U*par.Ai*(Tr-Tave)))/((mj*cpw)+(par.m_met*par.cp_met)));
F=dx1;
end
2 comentarios
bym
el 12 de Abr. de 2011
can you format all of the code using the {}code button
Andrew Newell
el 13 de Abr. de 2011
What is the "mistake", exactly? An error message?
Respuestas (0)
Categorías
Más información sobre Ordinary Differential Equations en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!