ode45, deval

Can't we know the function of yy values for successive xx values rather than the yy solution according to xx?

clear;

X0 = [0.1; 0.1; 0.1; 0.1];

no = 3712;

no_data = 3712;

H_total = 0;

time = load('time_epoxy.m');

dqdt = load('Q_epoxy.m');

temp = load('temp_epoxy.m');

baseline = zeros(no,1);

alpha = zeros(no,1);

dadt = zeros(no,1);

gradient = ( dqdt(no,1) - dqdt(1,1) ) / ( time(no,1) - time(1,1) );

for i = 1 : no

baseline(i,1) = dqdt(1,1) + gradient * ( time(i,1) - time(1,1) );

end

for i = 2 : no

H_total = H_total + ( dqdt(i,1) - baseline(i,1) ) * ( time(i,1) - time(i-1,1) ) ;

end

for i = 2 : no

alpha(i,1) = alpha(i-1,1) + ( ( dqdt(i,1) - baseline(i,1) ) * ( time(i,1) - time(i-1,1) ) ) / H_total;

end

for i = 1 : no

dqdt(i,1) = dqdt(i,1) - baseline(i,1);

end

for i = 1 : no

dadt(i,1) = dqdt(i,1) / H_total;

end

for i = 1 : no_data

X(i,1) = alpha( (i-1)*1+1 ,1);

Y(i,1) = dadt( (i-1)*1+1 ,1);

end

X(1,1) = alpha(1,1)+0.00001;

X(no_data,1) = alpha(no,1)-0.00001;

Y(no_data,1) = dadt(no,1);

options = optimset ('Largescale','off');

x=lsqnonlin(@fitting_epoxy_isothermal,X0,[],[],options,X,Y);

xt=real(x);

x=xt;

Y_sim = ( x(1) + x(2) * (alpha).^ x(3)).* ( 1-alpha).^ (x(4));

alpha_sim(1,1) = 0;

for i = 2 : no

alpha_sim(i,1) = alpha_sim(i-1,1) + Y_sim(i,1) * ( time(i,1) - time(i-1,1) ) ;

end

% 두번째 heating rate fitting

X0_2 = [0.1; 0.1; 0.1; 0.1];

no_2 = 2157;

no_data_2 = 2157;

H_total_2 = 0;

time_2 = load('time_epoxy_2.m');

dqdt_2 = load('Q_epoxy_2.m');

temp_2 = load('temp_epoxy_2.m');

baseline_2 = zeros(no_2,1);

alpha_2 = zeros(no_2,1);

dadt_2 = zeros(no_2,1);

gradient_2 = ( dqdt_2(no_2,1) - dqdt_2(1,1) ) / ( time_2(no_2,1) - time_2(1,1) );

for i = 1 : no_2

baseline_2(i,1) = dqdt_2(1,1) + gradient_2 * ( time_2(i,1) - time_2(1,1) );

end

for i = 2 : no_2

H_total_2 = H_total_2 + ( dqdt_2(i,1) - baseline_2(i,1) ) * ( time_2(i,1) - time_2(i-1,1) ) ;

end

for i = 2 : no_2

alpha_2(i,1) = alpha_2(i-1,1) + ( ( dqdt_2(i,1) - baseline_2(i,1) ) * ( time_2(i,1) - time_2(i-1,1) ) ) / H_total_2;

end

for i = 1 : no_2

dqdt_2(i,1) = dqdt_2(i,1) - baseline_2(i,1);

end

for i = 1 : no_2

dadt_2(i,1) = dqdt_2(i,1) / H_total_2;

end

for i = 1 : no_data_2

X_2(i,1) = alpha_2( (i-1)*1+1 ,1);

Y_2(i,1) = dadt_2( (i-1)*1+1 ,1);

end

X_2(1,1) = alpha_2(1,1)+0.00001;

X_2(no_data_2,1) = alpha_2(no_2,1)-0.00001;

Y_2(no_data_2,1) = dadt_2(no_2,1);

options = optimset ('Largescale','off');

x_2=lsqnonlin(@fitting_epoxy_isothermal,X0_2,[],[],options,X_2,Y_2);

xt_2=real(x_2);

x_2=xt_2;

Y_sim_2 = ( x_2(1) + x_2(2) * (alpha_2).^ x_2(3)).* ( 1-alpha_2).^ (x_2(4));

alpha_sim_2(1,1) = 0;

for i = 2 : no_2

alpha_sim_2(i,1) = alpha_sim_2(i-1,1) + Y_sim_2(i,1) * ( time_2(i,1) - time_2(i-1,1) ) ;

end

% 세번째 heating rate fitting

X0_3 = [0.1; 0.1; 0.1; 0.1];

no_3 = 907;

no_data_3 = 907;

H_total_3 = 0;

time_3 = load('time_epoxy_3.m');

dqdt_3 = load('Q_epoxy_3.m');

temp_3 = load('temp_epoxy_3.m');

baseline_3 = zeros(no_3,1);

alpha_3 = zeros(no_3,1);

dadt_3 = zeros(no_3,1);

gradient_3 = ( dqdt_3(no_3,1) - dqdt_3(1,1) ) / ( time_3(no_3,1) - time_3(1,1) );

for i = 1 : no_3

baseline_3(i,1) = dqdt_3(1,1) + gradient_3 * ( time_3(i,1) - time_3(1,1) );

end

for i = 2 : no_3

H_total_3 = H_total_3 + ( dqdt_3(i,1) - baseline_3(i,1) ) * ( time_3(i,1) - time_3(i-1,1) ) ;

end

for i = 2 : no_3

alpha_3(i,1) = alpha_3(i-1,1) + ( ( dqdt_3(i,1) - baseline_3(i,1) ) * ( time_3(i,1) - time_3(i-1,1) ) ) / H_total_3;

end

for i = 1 : no_3

dqdt_3(i,1) = dqdt_3(i,1) - baseline_3(i,1);

end

for i = 1 : no_3

dadt_3(i,1) = dqdt_3(i,1) / H_total_3;

end

for i = 1 : no_data_3

X_3(i,1) = alpha_3( (i-1)*1+1 ,1);

Y_3(i,1) = dadt_3( (i-1)*1+1 ,1);

end

X_3(1,1) = alpha_3(1,1)+0.00001;

X_3(no_data_3,1) = alpha_3(no_3,1)-0.00001;

Y_3(no_data_3,1) = dadt_3(no_3,1);

options = optimset ('Largescale','off');

x_3=lsqnonlin(@fitting_epoxy_isothermal,X0_3,[],[],options,X_3,Y_3);

xt_3=real(x_3);

x_3=xt_3;

Y_sim_3 = ( x_3(1) + x_3(2) * (alpha_3).^ x_3(3)).* ( 1-alpha_3).^ (x_3(4));

alpha_sim_3(1,1) = 0;

for i = 2 : no_3

alpha_sim_3(i,1) = alpha_sim_3(i-1,1) + Y_sim_3(i,1) * ( time_3(i,1) - time_3(i-1,1) ) ;

end

temp_peak=zeros(3,1); %% 각 heating rate 에서의 Tp

for i = 2:no

if dadt(i,1)>dadt(i-1,1)

temp_peak(1,1)=temp(i,1);

end

end

for i=2:no_2

if dadt_2(i,1)>dadt_2(i-1,1)

temp_peak(2,1)=temp_2(i,1);

end

end

for i=2:no_3

if dadt_3(i,1)>dadt_3(i-1,1)

temp_peak(3,1)=temp_3(i,1);

end

end

k1=[x(1); x_2(1); x_3(1)]; % Ea fitting 을 위한, 각 heating rate 에서의 k1

k2=[x(2); x_2(2); x_3(2)];

p1=polyfit(1000./(temp_peak+273.15),log(k1),1);

p2=polyfit(1000./(temp_peak+273.15),log(k2),1);

E1=p1(1); % linear regression 기울기 (-Ea/R)

E2=p2(1);

A1=exp(p1(2));

A2=exp(p2(2));

% E1=-E1*8.314472; % -E/R * 1000/Tp 에서 E 분리 (kJ/mol)

% E2=-E2*8.314472;

% E1 % Ea (kJ/mol)

c_x=1000./(temp_peak+273.15); % 그냥 curve fitting 앱 확인용 변수

c_y=log(k1);

sol = ode45(@(t,y) x(1)+x(2)*y^x(3)*(1-y)^x(4), [12 24], 0); % da/dt -> a 수치적분

xx = linspace(12,24,10000);

yy = deval(sol,xx);

figure(4)

plot(xx,yy)

grid on;