Integrator_test.m

clear;
clc;
close all;
hold off;

dt = 1.19;
time = 50;
n_steps = round(time/dt);


m = 500;
c = 0.5;
h = 410;

T0 = 475;
Q0 = m*c*T0;

T_ext = 280;

Texp = [];
Qexp = [];
Qdotexp = [];
texp = [];

Texp(1) = T0;
Qexp(1) = Q0;
Qdotexp(1) = 0;
texp(1) = 0;

for ii = 1:n_steps
  Qdotexp(ii) = h*(T_ext-Texp(ii));
  Qexp(ii+1) = Qexp(ii) + Qdotexp(ii)*dt;
  Texp(ii+1) = Qexp(ii+1)/(m*c);
  texp(ii+1) = texp(ii)+dt;
endfor

plot(texp,Texp)
hold on;


Timp(1) = T0;
Qimp(1) = Q0;
Qdotimp(1) = 0;
timp(1) = 0;

for ii = 1:n_steps
  
  Qimp(ii+1) = Qimp(ii);
  
  Qimp(ii+1) = Qimp(ii) + h*(T_ext-(Qimp(ii+1)/(m*c)))*dt;
  
  for jj = 1:10
    Qimp(ii+1) = Qimp(ii) + h/2*(((T_ext-(Qimp(ii+1)/(m*c))))+((T_ext-(Qimp(ii)/(m*c)))))*dt;
    
  endfor

  Timp(ii+1) = Qimp(ii+1)/(m*c);
  
  timp(ii+1) = timp(ii)+dt;
endfor
plot(timp,Timp)


Timp(1) = T0;
Qimp(1) = Q0;
Qdotimp(1) = 0;
timp(1) = 0;

for ii = 1:n_steps
  
  Qimp(ii+1) = Qimp(ii);
  
  Qimp(ii+1) = Qimp(ii) + h*(T_ext-(Qimp(ii+1)/(m*c)))*dt;
  
  for jj = 1:50
    Qimp(ii+1) = Qimp(ii) + h/2*(((T_ext-(Qimp(ii+1)/(m*c))))+((T_ext-(Qimp(ii)/(m*c)))))*dt;
    
  endfor

  Timp(ii+1) = Qimp(ii+1)/(m*c);
  
  timp(ii+1) = timp(ii)+dt;
endfor
plot(timp,Timp)

Qeq = [];
Qdoteq = [];
Teq = [];
teq = [];

Teq(1) = T0;
teq(1) = 0;
Qdoteq(1) = 0;
Qeq(1) = Q0;

for ii = 1:n_steps
  Qdoteq(ii+1) = m*c*(Teq(ii)-T_ext)*(1-exp(-(h*dt)/(m*c)));
  Qeq(ii+1) = Qeq(ii)-Qdoteq(ii+1);
  
  Teq(ii+1) = Qeq(ii+1)/(m*c);
  teq(ii+1) = teq(ii)+dt;
endfor
plot(teq,Teq)

hold off;
legend("Euler","Implicit Trapazoid n = 10",...
 "Implicit Trapazoid n = 50","Piecewise Analytical");
 
xlabel("Time [sec]");
ylabel("Temperature [K]");
title("Integrator Test")