ODE 45 and function combination

For my study I recieved an assignment to calculate the most optimal way of isolating a house. This is done using matlab.

The assignment is included underneath.

Me and my group are getting several errors, and really dont know what is going wrong.

We would really appreciate it if someone could help us out.

The script

% define constants
A = 1;
B = 2;
Cint = 1.0*10^4;
Cwall = 5.0*10^2;
R1 = 0.1;
R2 = 0.08;
Rwin = 0.1;
Tamb = dlmread('Temperaturesummer30.dat');
QHVAC = 0;

% time dependent window 
tspan = [0:1:1440];  % calculate from t=0 up to t=1440
                % can define the time interval steps
y0 = [18 20];    % initial conditions
Twall = y0(1);
Tint = y0(2);


% define your ode-function
odefunc = @(t,dydt) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

% integrate the system of differential equations from tspan(1) to 
% tspan(2) with initial conditions y0

[t,dydt] = ode45(odefunc, tspan, y0);

% plot the results
plot(t,dydt)

The function

function [dydt] = odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC)
% set of ordinary differential equations
% input:  A - constant
%         B - constant
%         t - the time variable
%         y - the state variable
%output:  dydt - the set of equations



% initialize the set of equations
dydt = zeros(2,1);

% define the set of equations
dydt(1) = ((Tamb-Twall)/(Cwall *R2)) + ((Tint-Twall)/(Cwall*R1)) %Twall                       %y(2)+A;
dydt(2) = ((Twall-Tint))/(Cint *R1) + ((Tamb-Tint)/(Cint* Rwin)) + (QHVAC/Cint) %Tint        %(A/B)*y(1);
end



The errors 

Unable to perform assignment because the left and right sides have a different number of elements.

Error in odeFunction (line 15)
dydt(1) = ((Tamb-Twall)/(Cwall *R2)) + ((Tint-Twall)/(Cwall*R1)) %Twall                       %y(2)+A;

Error in odeScript>@(t,dydt)odeFunction(t,Twall,Tint,Tamb,Cwall,Cint,R1,R2,Rwin,QHVAC) (line 26)
odefunc = @(t,dydt) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in odeScript (line 31)
[t,dydt] = ode45(odefunc, tspan, y0);
 
>> 

[Ameer Hamza]

Where is y vector inside odeFunction?

[Matthijs neerhof]

Thanks for your answer!

I had this at first, but it still did not work. To experiment I removed this. I will add it again.

[Matthijs neerhof]

I assume you mean y0 between the brackets at the top line of the function?

[Ameer Hamza]

No, I mean that inside the following function, you haven't used the state variable 'y' or time 't' anywhere

function [dydt] = odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC)
% set of ordinary differential equations
% input:  A - constant
%         B - constant
%         t - the time variable
%         y - the state variable
%output:  dydt - the set of equations
% initialize the set of equations
dydt = zeros(2,1);
% define the set of equations
dydt(1) = ((Tamb-Twall)/(Cwall *R2)) + ((Tint-Twall)/(Cwall*R1)) %Twall                       %y(2)+A;
dydt(2) = ((Twall-Tint))/(Cint *R1) + ((Tamb-Tint)/(Cint* Rwin)) + (QHVAC/Cint) %Tint        %(A/B)*y(1);
end

Considering above, the line

odefunc = @(t,dydt) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

seems incorrect. Can you show your ODE in mathematical form?

[Matthijs neerhof]

What do you mean with mathematical form?

Answer 1

[Star Strider]

I cannot run your code, since I do not have the ‘Tamb’ data.

Nevertheless ‘Tamb’ is likely the problem, since it appears to be an array of some sort, and that is throwing the error. If it is time-dependent, it will be necessary to interpolate it using the approach in ODE with Time-Dependent Terms in the ode45 (and related functions) documentation.

[Matthijs neerhof]

Thanks!

I now added Temp = interp1(1:1441, Tamb, t) in the function,

and I changed Tamb in the equations in the function to Temp.

I am now getting results without errors. Nevertheless, the graph that pops up is empty.

Any idea on how to solve this?

[Star Strider]

It is a bit difficult to determine what the problem is when I cannot run your code.

I suspect ‘odefunc’ should be:

odefunc = @(t,Twall) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

since I assume ‘Twall’ is what you want, however it is a bit difficult to follow your code.

Note that ‘dydt’ is not an argument to either ‘odefunc’ or ‘odeFunction’, so I am a bit surprised that did not itself throw an error.

[Matthijs neerhof]

The result should be a graph with Tamb, Twall, and Tint.

temperature on the y-axis and time on the x-axis.

Right now the graph seems to be empty. I removed dydt from both places.

I would like to send you the temperaturesummer30 file, but it is a .dat file which is not supported.

[Star Strider]

The best way to send it would be to read it and save it (or at least all the relevant variables) to a .mat file, then attach that.

[Matthijs neerhof]

Thanks!

I attatched the file, I hope it works!

[Star Strider]

My pleasure!

All the variables appear to be there, so I will spend some time with it in a few minutes.

You mentioned that you updated your code, so I would appreciate you posting the updated code as well.

One problem is that when I load it, this Warning appears:

Warning: Could not find appropriate function on path loading function handle X:\My
Downloads\odeScript.m>@(t,dydt)odeFunction(t,Twall,Tint,Tamb,Cwall,Cint,R1,R2,Rwin,QHVAC) 

No worries since all the variables appear to be present.

[Matthijs neerhof]

script
% script for solving an ode45
clear
clc

% define constants
A = 1;
B = 2;
Cint = 1.0*10^4;
Cwall = 5.0*10^2;
R1 = 0.1;
R2 = 0.08;
Rwin = 0.1;
Tamb = dlmread('Temperaturesummer30.dat');
QHVAC = 0;

% time dependent window 
tspan = [0:1:1440];  % calculate from t=0 up to t=1440
                % can define the time interval steps
y0 = [18 20];    % initial conditions
Twall = y0(1);
Tint = y0(2);


% define your ode-function
odefunc = @(t,dydt) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

% integrate the system of differential equations from tspan(1) to 
% tspan(2) with initial conditions y0

[t,dydt] = ode45(odefunc, tspan, y0);

% plot the results
plot(t,dydt)



%function
function [dydt] = odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC)
% set of ordinary differential equations
% input:  A - constant
%         B - constant
%         t - the time variable
%         y - the state variable
%output:  dydt - the set of equations



% initialize the set of equations
dydt = zeros(2,1);

% define the set of equations
Temp = interp1(1:1441, Tamb, t) ;
dydt(1) = ((Temp-Twall)/(Cwall *R2)) + ((Tint-Twall)/(Cwall*R1)) %Twall                       %y(2)+A;
dydt(2) = ((Twall-Tint))/(Cint *R1) + ((Temp-Tint)/(Cint* Rwin)) + (QHVAC/Cint) %Tint        %(A/B)*y(1);
end




%thanks again for your help!

[Star Strider]

My pleasure!

One problem is that it is necessary to reference the appropriate elements of ‘Twall’ in the differential equation. I used ‘Twall(1)’ in ‘dydt(1)’ and ‘Twall(2)’ in ‘dydt(2)’. I have no idea if that is correct, so you will need to sort that and provide the correct references in each equation. Nevertheless, with these changes, your code runs without error.

The only other change I made in ‘odeFunction’ is:

Temp = interp1(0:1440, Tamb, t) ;

(note that it uses (0:1440) instead of (1:1441)).

The lightly-edited function is now:

function [dydt] = odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC)
% set of ordinary differential equations
% input:  A - constant
%         B - constant
%         t - the time variable
%         y - the state variable
%output:  dydt - the set of equations



% initialize the set of equations
dydt = zeros(2,1);

% define the set of equations
Temp = interp1(0:1440, Tamb, t) ;
dydt(1) = ((Temp-Twall(1))/(Cwall *R2)) + ((Tint-Twall(1))/(Cwall*R1)) %Twall                       %y(2)+A;
dydt(2) = ((Twall(2)-Tint))/(Cint *R1) + ((Temp-Tint)/(Cint* Rwin)) + (QHVAC/Cint) %Tint        %(A/B)*y(1);
end

Changes in your script:

% define your ode-function
odefunc = @(t,Twall) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

% integrate the system of differential equations from tspan(1) to 
% tspan(2) with initial conditions y0

[t,y] = ode45(odefunc, tspan, y0);

% plot the results
figure
plot(t,y)
grid

This now runs without error.

I defer to you to be certain the coding logic is correct in order to give the correct results.

[Matthijs neerhof]

Thanks!

I am again getting a graph with no results. The results should be as the file I added. I hope you can help us one last time.

[Matthijs neerhof]

I now also found a way to share the .dat file.

If I am correct it is in the zip-file I attatched.

[Star Strider]

I get this plot:

It is not exactly the one in the assignment, however you will need to sort that, since it most likely depends on where the 2 values of ‘Twall’ belong, since I am not certain where ‘Twall(1)’ and ‘Twall(2)’ go in the derivatives. I assigned them as I did simply to get the code to run without throwing errors.

The ‘odeFunction’ file did not change from my previous Comment. The calls to it and the plot call are:

% define your ode-function
odefunc = @(t,Twall) odeFunction(t, Twall, Tint, Tamb, Cwall, Cint, R1, R2, Rwin, QHVAC);

% integrate the system of differential equations from tspan(1) to 
% tspan(2) with initial conditions y0

[t,y] = ode45(odefunc, tspan, y0);

% plot the results
figure
plot(t,y)
grid
xlabel('time (in minutes)')
ylabel('Temperature (celsius)')

as previously (with only the axis titles added).