ode45, deval

1 次查看（过去 30 天）
Yongchul Shin
Yongchul Shin 2019-9-27
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;

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by