Defining Differential Equation in Runge Kutta 4th order

Hello!

I am trying to solve the following differential equation with runge-kutta, 4th order but I have no idea how to define the equation in matlab :/

What I have so far is:

%functional terms

Pel= rho*l*I^(2)/Aw;

Pc = h*Al*(T-Ta);

Pe = boltzmann*e*al*(T^(4)-Ta^(4));

%Differential equation

h=0.1; a=0; b=100; % h is the step size, t=[a,b] t-range
t = a:h:b; % Computes t-array up to t=100
T = zeros(1,numel(t)); % Memory preallocation
T(0) = Ta; % initial condition; in MATLAB indices start at Ta
T_dot =@(T)(delta/(m*ct))*Pel(T)-(Pc(T)+Pe(T)); %insert function to be solved
% the function is the expression after (T,t)
% table title
fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','t(i)','k1','k2','k3',
'k4','y(i)');
for ii=1:1:numel(t)
 k1 = T_dot(t(ii),T(ii));
 k2 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k1);
 k3 = T_dot((t(ii)+0.5*h),(T(ii)+0.5*h*k2));
 k4 = T_dot(t(ii)+h),(T(ii)+h*k3));
 T(ii+1) = T(ii) + (h/6)*(k1+2*k2+2*k3+k4); % main equation
 % table data
 fprintf(fid,'%7d %7.2f %7.3f %7.3f',ii, t(ii), k1, k2);
 fprintf(fid,' %7.3f %7.3f %7.3f \n', k3, k4, T(ii));
end
T(numel(t))=[ ]; % erase the last computation of T(n+1)
% Solution PLOT:
plot(t,T,' ok')
title('RK-4--Numerical Solution---');
ylabel('T'); xlabel('t'); legend('numerical');
grid on
fclose(fid);

Honestly, I have copied the code from a site and tried to adapt it to my equation, but as I am a complete beginner I have no idea where or what I did wrong, I just know that it's wrong.

Maybe someone can help me out? I'd appreciate every help :)!

[Torsten]

We need rho, l, I, Aw, h, Ta, boltzmann, e, al, delta, m, ct.

And your function definition is not the same as in your mathematical formula.

[phkstudent]

I defined those parameters before the given code, but I left it out cause I thought it might look cinfusing and too much

and yes, I have noticed that I solved my equation incorrectly (though Pw is left out on purpose)

Ta = 300; 
g = 9.81;
V=12; %Volt
I=10; %Ampere

delta = 1; %switch
r=0.01; %[m]
l=0.2; %[m]
Aw = pi*r^2;
Al = l*2*r;
m=6500*Aw*l;
e=0.37;
ct= 470;
boltzmann = 5.67*10^(-8); %[W/m^2K^4]

kinvis= 15.89;
dynvis= 184.6*10^(-13);
cp= 1.007; %specific heat
kg= 26.3^(-6); %thermal conductivity gas
Pr = dynvis*cp/kg;

Gr = (2*(T-Ta)/(T+Ta))*g*r^(3)/(8*kinvis^(2));

Nu = (1+0.287*((Gr*Pr)/(1+(0.56/Pr)^(9/16))^(9/16)))^(2); %Nusselt

Kg = 0.026;

rho = V*Aw/(I*l);
h=Nu*kg/l; %thermal conductivity coefficent

Answer 1

[Torsten]

Ta = 300; 
g = 9.81;
V=12; %Volt
I=10; %Ampere
delta = 1; %switch
r=0.01; %[m]
l=0.2; %[m]
Aw = pi*r^2;
Al = l*2*r;
m=6500*Aw*l;
e=0.37;
ct= 470;
boltzmann = 5.67*10^(-8); %[W/m^2K^4]
kinvis= 15.89;
dynvis= 184.6*10^(-13);
cp= 1.007; %specific heat
kg= 26.3^(-6); %thermal conductivity gas
Pr = dynvis*cp/kg;
Gr = @(T)(2*(T-Ta)/(T+Ta))*g*r^(3)/(8*kinvis^(2));
Nu = @(T)(1+0.287*((Gr(T)*Pr)/(1+(0.56/Pr)^(9/16))^(9/16)))^(2); %Nusselt
Kg = 0.026;
rho = V*Aw/(I*l);
h=@(T)Nu(T)*kg/l; %thermal conductivity coefficent
%functional terms
Pel= @(T)rho*l*I^(2)/Aw;
Pc = @(T)h(T)*Al*(T-Ta);
Pe = @(T)boltzmann*e*Al*(T^4-Ta^4);
Pw = @(T)0;
%Differential equation
h=0.1; a=0; b=100; % h is the step size, t=[a,b] t-range
t = a:h:b; % Computes t-array up to t=100
T = zeros(1,numel(t)); % Memory preallocation
T(1) = Ta; % initial condition; in MATLAB indices start at Ta
T_dot =@(t,T)(delta*Pel(T)-(Pc(T)+Pe(T)+delta*Pw(T)))/(m*ct); %insert function to be solved
% the function is the expression after (T,t)
% table title
%fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','t(i)','k1','k2','k3',
%'k4','y(i)');
for ii=1:1:numel(t)-1
    k1 = T_dot(t(ii),T(ii));
    k2 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k1);
    k3 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k2);
    k4 = T_dot(t(ii)+h,T(ii)+h*k3);
    T(ii+1) = T(ii) + (h/6)*(k1+2*k2+2*k3+k4); % main equation
    % table data
    %fprintf(fid,'%7d %7.2f %7.3f %7.3f',ii, t(ii), k1, k2);
    %fprintf(fid,' %7.3f %7.3f %7.3f \n', k3, k4, T(ii));
end
%T(numel(t))=[ ]; % erase the last computation of T(n+1)
% Solution PLOT:
plot(t,T,' ok')
title('RK-4--Numerical Solution---');
ylabel('T'); xlabel('t'); legend('numerical');
grid on

[phkstudent]

WOW! Thank you so much!! <3

[phkstudent]

Hey, I do have another question if thats okay :D

If the current I would also be a variable, meaning that T = f(t,I), how can I plot this?

I wanted I to be like an axis as well, but I am not even sure if RK4 method is even possible when I have 2 depending variables.

[Torsten]

Do you want to solve the equation for different (constant) values for I (thus get several (t,T) solution curves) or do you want to change I over time in one calculation (thus get one (t,T) solution curve with time-dependent I) ?

[phkstudent]

I want to change it to T over I, and if possible (t,T) curve with depency on I

[Torsten]

So in each simulation, "I" can remain constant over time, but you want to solve the equation for different values for "I" ?

[phkstudent]

By the end I want to display which current "I" do I need to reach a specific temperature T and in what time is that temperature reached

[Torsten]

%Differential equation
h=0.1; a=0; b=100; % h is the step size, t=[a,b] t-range
t = a:h:b; % Computes t-array up to t=100
T = zeros(1,numel(t)); % Memory preallocation
I_array = 5:5:100;
T_array = zeros(numel(I_array),numel(t));
for i = 1:numel(I_array)
    Ta = 300; 
    g = 9.81;
    V=12; %Volt
    %I=10; %Ampere
    I = I_array(i); % Ampere
    delta = 1; %switch
    r=0.01; %[m]
    l=0.2; %[m]
    Aw = pi*r^2;
    Al = l*2*r;
    m=6500*Aw*l;
    e=0.37;
    ct= 470;
    boltzmann = 5.67*10^(-8); %[W/m^2K^4]
    kinvis= 15.89;
    dynvis= 184.6*10^(-13);
    cp= 1.007; %specific heat
    kg= 26.3^(-6); %thermal conductivity gas
    Pr = dynvis*cp/kg;
    Gr = @(T)(2*(T-Ta)/(T+Ta))*g*r^(3)/(8*kinvis^2);
    Nu = @(T)(1+0.287*((Gr(T)*Pr)/(1+(0.56/Pr)^(9/16))^(9/16)))^2; %Nusselt
    Kg = 0.026;
    rho = V*Aw/(I*l);
    h_value=@(T)Nu(T)*kg/l; %thermal conductivity coefficent
    %functional terms
    Pel= @(T)rho*l*I^2/Aw;
    Pc = @(T)h_value(T)*Al*(T-Ta);
    Pe = @(T)boltzmann*e*Al*(T^4-Ta^4);
    Pw = @(T)0;
    %Differential equation
    %h=0.1; a=0; b=100; % h is the step size, t=[a,b] t-range
    %t = a:h:b; % Computes t-array up to t=100
    %T = zeros(1,numel(t)); % Memory preallocation
    T(1) = Ta; % initial condition; in MATLAB indices start at Ta
    T_dot =@(t,T)(delta*Pel(T)-(Pc(T)+Pe(T)+delta*Pw(T)))/(m*ct); %insert function to be solved
    % the function is the expression after (T,t)
    % table title
    %fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','t(i)','k1','k2','k3',
    %'k4','y(i)');
    for ii=1:numel(t)-1
        k1 = T_dot(t(ii),T(ii));
        k2 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k1);
        k3 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k2);
        k4 = T_dot(t(ii)+h,T(ii)+h*k3);
        T(ii+1) = T(ii) + (h/6)*(k1+2*k2+2*k3+k4); % main equation
        % table data
        %fprintf(fid,'%7d %7.2f %7.3f %7.3f',ii, t(ii), k1, k2);
        %fprintf(fid,' %7.3f %7.3f %7.3f \n', k3, k4, T(ii));
    end
    T_array(i,:) = T;
end
%T(numel(t))=[ ]; % erase the last computation of T(n+1)
% Solution PLOT:
plot(t,T_array,' ok')
title('RK-4--Numerical Solution---');
ylabel('T'); xlabel('t'); legend('numerical');
grid on

[phkstudent]

oh okay, so each curve represents the T-t solution for a specific curren?

Thank you so much!!

[Torsten]

Or even better:

%Differential equation
h=0.1; a=0; b=100; % h is the step size, t=[a,b] t-range
t = a:h:b; % Computes t-array up to t=100
I = 5:5:100;
T = zeros(numel(I),numel(t));
for i = 1:numel(I)
    T(i,:) = Temperature(t,I(i));
end
plot(t,T,' ok')
title('RK-4--Numerical Solution---');
ylabel('T'); xlabel('t'); legend('numerical');
grid on

function [T] = Temperature(t,I)
  Ta = 300; 
  g = 9.81;
  V=12; %Volt
  %I=10; %Ampere
  delta = 1; %switch
  r=0.01; %[m]
  l=0.2; %[m]
  Aw = pi*r^2;
  Al = l*2*r;
  m=6500*Aw*l;
  e=0.37;
  ct= 470;
  boltzmann = 5.67*10^(-8); %[W/m^2K^4]
  kinvis= 15.89;
  dynvis= 184.6*10^(-13);
  cp= 1.007; %specific heat
  kg= 26.3^(-6); %thermal conductivity gas
  Pr = dynvis*cp/kg;
  Gr = @(T)(2*(T-Ta)/(T+Ta))*g*r^(3)/(8*kinvis^(2));
  Nu = @(T)(1+0.287*((Gr(T)*Pr)/(1+(0.56/Pr)^(9/16))^(9/16)))^(2); %Nusselt
  Kg = 0.026;
  rho = V*Aw/(I*l);
  h_value=@(T)Nu(T)*kg/l; %thermal conductivity coefficent
  %functional terms
  Pel= @(T)rho*l*I^(2)/Aw;
  Pc = @(T)h_value(T)*Al*(T-Ta);
  Pe = @(T)boltzmann*e*Al*(T^4-Ta^4);
  Pw = @(T)0;
  %Differential equation
  %h=0.1; a=0; b=100; % h is the step size, t=[a,b] t-range
  %t = a:h:b; % Computes t-array up to t=100
  T = zeros(1,numel(t)); % Memory preallocation
  T(1) = Ta; % initial condition; in MATLAB indices start at Ta
  T_dot =@(t,T)(delta*Pel(T)-(Pc(T)+Pe(T)+delta*Pw(T)))/(m*ct); %insert function to be solved
  % the function is the expression after (T,t)
  % table title
  %fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','t(i)','k1','k2','k3',
  %'k4','y(i)');
  for ii=1:1:numel(t)-1
    h = t(ii+1)-t(ii);
    k1 = T_dot(t(ii),T(ii));
    k2 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k1);
    k3 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k2);
    k4 = T_dot(t(ii)+h,T(ii)+h*k3);
    T(ii+1) = T(ii) + (h/6)*(k1+2*k2+2*k3+k4); % main equation
    % table data
    %fprintf(fid,'%7d %7.2f %7.3f %7.3f',ii, t(ii), k1, k2);
    %fprintf(fid,' %7.3f %7.3f %7.3f \n', k3, k4, T(ii));
  end
  %T(numel(t))=[ ]; % erase the last computation of T(n+1)
  % Solution PLOT:
  %plot(t,T,' ok')
  %title('RK-4--Numerical Solution---');
  %ylabel('T'); xlabel('t'); legend('numerical');
  %grid on
end

[phkstudent]

This looks very good, thanks!

[phkstudent]

I have changed some constant values into values depending on T, but now it gives me this error:

Unable to perform assignment because the size of the left side is 1-by-201 and the size of the right side

is 1-by-51.

Error in tempbehaviorVScurrent (line 17)

T(i,:) = Temperature(t,I(i));

%Differential equation
h=0.5; a=0; b=100; % h is the step size, t=[a,b] t-range
t = a:h:b; % Computes t-array up to t=100
I = 3:1:10;
T = zeros(numel(I),numel(t));
for i = 1:numel(I)
    T(i,:) = Temperature(t,I(i));
end

plot(t,T,' ok')
title('RK-4--Numerical Solution---');
ylabel('T'); xlabel('t'); legend('numerical');
grid on

function [T] = Temperature(t,I)
Ta = 296.15;                     %ambient Temp in [K]
g = 9.81;                        %[m^2/s]
V=12;                            %voltage [V]                             
delta = 1;                       %switch
r=0.0015;                        %[m]
l= 0.1;                          %[m]
Aw = pi*r^2;                     %[m^2]
Al = l*2*r;                      %Schnittfläche längs[m^2]
Au = pi*2*r*l/2;                 %1/2 Umfangsfläche längs [m^2]
density_w = 6500;                %[kg/m^3]
m=density_w*Aw*l;                %[kg]
rho = V*Aw/(I*l);                %wire resistivity
ew=0.37;                         %emissivity of wire
e_mdf=0.85;                      %emissivity of MDF
ct= 470;                         %total heat capacity [J/kgK]
boltzmann = 5.67*10^(-8);        %Boltzmann constant [W/m^2K^4]
p = 1.013*10^5;                  %air pressure [Pa]
R = 0.287/1000;                  %gas constant [J/kgK]
Z = 1;                           %compressibility factor, ideal gas

%Terms depending on T
density_air = @(T) (p/(R*((T+Ta)/2)*Z));                  %[kg/m^2]
dynvis= @(T) (0.9228*(5*10^(-8)*((T+Ta)/2)+5*10^(-6)));   %[kg/ms]
kinvis= @(T) dynvis(T)/density_air(T);                    %kinematic viscosity air [m^2/s)

%Gr = @(T)(2*(T-Ta)/(T+Ta))*g*r^(3)/(8*kinvis(T)^(2));                  &Grashof 
%Pr = @(T) (dynvis(T)*cp/kg(T));                                        %Prandtl
%Nu = @(T)(1+0.287*((Gr(T)*Pr(T))/(1+(0.56/Pr(T))^(9/16))^(9/16)))^(2); %Nusselt

cp= 1003;                                                %specific heat of air [J/kgK]
kg= @(T) (1.04*(7*10^(-5)*((T+Ta)/2)+0.0037));            %gas thermal conductivity [W/mK]
lambda = 0.0262;                                          %wärmeleitfähigkeit air [W/mK]

Pr = @(T) (dynvis(T)*cp/lambda);
Gr = @(T) g*(T-Ta)*l^(3)/(Ta*kinvis(T)^(2));
Ra = @(T) Gr(T)*Pr(T);
Nu = @(T) 1.158*Ra(T)^(0.1196);
h_value=@(T)Nu(T)*kg(T)/l;                             %thermal conductivity coefficent


               
%functional terms of differential equation
Pel= @(T)rho*l*I^(2)/Aw;                      %input power
QconvA = @(T)h_value(T)*Au*(T-Ta);                  %heat convection wire->air
QradMDF = @(T) e_mdf*boltzmann*Au*(T^4-Ta^4); %radiation heat transfer wire->MDF
QradA = @(T)ew*boltzmann*Au*(T^4-Ta^4);       %radiation heat transfer wire->air


%Runge-Kutta 4th order
%h=1; a=0; b=50;                    % h is the step size, t=[a,b] time-range in [s]
%t = a:h:b;                         
T = zeros(1,numel(t));             % Memory preallocation
T(1) = Ta;                         % initial condition; in MATLAB indices start at Ta
T_dot =@(t,T)(delta*Pel(T)-(QconvA(T)+QradA(T)+QradMDF(T)))/(m*ct); % differential equation to be solved

%fprintf(fid,'%7s %7s %7s %7s %7s %7s %7s \n','i','t(i)','k1','k2','k3','k4','y(i)');

for ii=1:1:numel(t)-1
    h = t(ii+1)-t(ii);
    k1 = T_dot(t(ii),T(ii));
    k2 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k1);
    k3 = T_dot(t(ii)+0.5*h,T(ii)+0.5*h*k2);
    k4 = T_dot(t(ii)+h,T(ii)+h*k3);
    T(ii+1) = T(ii) + (h/6)*(k1+2*k2+2*k3+k4); % main equation
    % table data
    %fprintf(fid,'%7d %7.2f %7.3f %7.3f',ii, t(ii), k1, k2);
    %fprintf(fid,' %7.3f %7.3f %7.3f \n', k3, k4, T(ii));
  end
%   T(numel(t))=[ ]; % erase the last computation of T(n+1)
%   Solution PLOT:
%   plot(t,T,' ok')
%   title('RK-4--Numerical Solution---');
%   ylabel('T'); xlabel('t'); legend('numerical');
%   grid on
end

[Torsten]

The heat transfer coefficient h had to be changed in h_value because h is the stepsize of the time integration.

Further, the t-array is already defined in the calling program - so I commented it out in the function.

[phkstudent]

You truly are the best! <3