Using Ode45 to solve differential equation with time dependent variable

Hi my name is Hidde and i am a first year IEM student,

I got an assignment where a part of it is to solve 2 differential equations and plot them in a graph. I just keep getting error and i cannot get the code to work.

I have some experience with matlab but i never used Ode45 before. I've read numerous tutorials and watched a bunch of video's but i really need help. Tamb is time dependant and changes every second.

OdeScript

% script for solving an ode45
clear all
clc
% define constants
R1 = 0.4;
R2 = 0.3;
Rwin = 0.5;

Cint = 1.0 * 10^4;
Cwall = 5.0 * 10^2; 

Tamb = dlmread('Team_82.dat');


A = 1 / (Cwall*R2);
B = 1 / (Cwall*R1);
C = 1 / (Cint*R1);
D = 1 / (Cint*Rwin);

% time dependent window 
tspan = linspace(0,1440,1441);  % calculate from t=0 up to t=1440

y0 = [18 20];    % initial conditions

% define your ode-function
odefunc = @(t,dTempdt) odeFunction1(t,Twall, Tint, Tamb, A, B, C, D);



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

% plot the results
plot(t,dTempdt)

OdeFunction

function [dTempdt] = odeFunction1(t,Twall, Tint, Tamb, A, B, C, D)
% set of ordinary differential equations
% input:  A - constant
%         B - constant
%         C - constant
%         D - constant
%         t - the time variable
%         Tint - state variable
%         Twall - state variable
%output:  dTempdt - the set of equations

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

% define the set of equations
dTempdt(1) = A*Tamb - A*Twall + B*Tint - B*Twall;
dTempdt(2) = C*Twall - C*Tint + D*Tamb - D*Tint;
end

I am getting the folowing errors:

Undefined function or variable 'Twall'.

Error in odeScript1>@(t,dTempdt)odeFunction1(t,Twall,Tint,Tamb,A,B,C,D)

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 odeScript1 (line 31)
[t,dTempdt] = ode45(odefunc, tspan, y0);

Any help would be greatly appreciated!

[Stephan]

please attach the needed data

[Hidde Kemperink]

Hi Stephan,

I attached the data and the assignment pdf in data.zip

[Daniel]

Hi Hidde,

you're making a mistake when defining the ode-function. The function call to the pseudo-function

odefunc = @(t,dTempdt) odeFunction1(t,Twall, Tint, Tamb, A, B, C, D);

accepts two input arguments, as it should. Namely the time and the state vector: (t,y). It has to output dy/dt. The problem is you define this pseudo-function "odefunc" with two input parameters (t and dTempdt) but you have it call another function "odeFunction1" with other input parameters (t,Twall,Tint,Tamb,A,B,C,D), hence the error message. odeFunction1 does not know the values for all those parameters/constants, if you call it this way.

You can only have one state vector y, so your state variables should be grouped inside one vector. Then the ode-function accepts two inputs (time t, state vector y) and needs to calculate dy/dt. To do that you need to define the respective equations inside this ode-function. If these equations use constants or other parameters, you either define them inside the function or use global variables for example.

So the structure should look similar to this:

odefunc = @(t,y) odeFunction1(t,y) %output of odeFunction1 has to be dy/dt

function [dydt] = odeFunction1(t,y)
R1 = 0.4;
R2 = 0.3;
Rwin = 0.5;
Cint = 1.0 * 10^4;
Cwall = 5.0 * 10^2;

A = 1 / (Cwall*R2);
B = 1 / (Cwall*R1);
C = 1 / (Cint*R1);
D = 1 / (Cint*Rwin);

Twall = y(1,1); %these are your state variables inside the state vector y
Tint = y(2,1); % you have two states (Twall and Tint)

Tamb = ...?; %You need to work out the current value for Tamb...
% in every time step. Hint: you can use the input parameter t (current time)...
% for that. Consider the task in your assignment paper and check the Matlab...
% Documentation

dydt = zeros(2,1); %your output time derivative of the state vector
dydt(1,1) = A*Tamb - A*Twall + B*Tint - B*Twall;
dydt(2,1) = C*Twall - C*Tint + D*Tamb - D*Tint;
end

Finally you have to make sure that the initial conditions are defined with the same dimensions as your state vector in the ode-function. So if you're using dydt = zeros(2,1) it is a 2x1 vector:

y0 = [18; 20];    % initial conditions

[t_solution,y_solution] = ode45(odefunc, tspan, y0);

Hope this helps... ;)

Cheers,

Daniel

[Hidde Kemperink]

Hi Daniel,

First of all thank you so much for your help! i just cant figure out how to get Tamb to change when the value for t changes. I tried numerous tutorials but i just cannot figure it out. If you could help me with that i would be very thankfull!

this is what i tried:

importing data as:

Data = dlmread('Team_82.dat');

and defining Tamb as, t+1 since the index starts at 1 and t starts at 0:

Tamb = Data(t+1);

[Daniel]

Ok, in this case do the following:

• Declare Tamb global

global Tamb;
Tamb = ...; %import your data here
    

• Inside odeFunction1 use an interpolation function that returns the current Temp-value for a give input time using linear interpolation. (As Stephan used in his answer). You also have to define the points in time for the data you import. (tspan)

global Tamb; %this ensures that the function can access the data you read...
% into the Matlab base workspace.
tspan = linspace(0,length(Tamb),length(Tamb)+1);
Tamb_t = interp1(tspan,Tamb,t);

[Hidde Kemperink]

Well the script & function work when i just assign a set value to Tamb. I tried the code for importing the data but i get there errors. I made sure to change Tamb to Tamb_t in the defined DE's

Error using interp1>reshapeAndSortXandV (line 424)
X and V must be of the same length.

Error in interp1 (line 93)
    [X,V,orig_size_v] = reshapeAndSortXandV(varargin{1},varargin{2});

Error in odeFunction1 (line 16)
    Tamb_t = interp1(tspan,Tamb,t);

Error in odeScript1>@(t,y)odeFunction1(t,y)

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 odeScript1 (line 20)
[t_solution,y_solution] = ode45(odefunc, tspan, y0);

[Hidde Kemperink]

Okay, so i changed

tspan = linspace(0,length(Tamb),length(Tamb)+1);

to

tspan = linspace(0,length(Tamb)-1,length(Tamb));

and now it plot a proper graph, im not sure what im doing is correct tho

[Stephan]

I changed the name - now it should be clear that i use ode45...

Answer 1

[Stephan]

% call the nested function
[t,dTempdt] = my_assignment;
% plot the results
plot(t,dTempdt)

function [t,dTempdt] = my_assignment
% define constants
R1 = 0.4;
R2 = 0.3;
Rwin = 0.5;
Cint = 1.0 * 10^4;
Cwall = 5.0 * 10^2; 
Tamb = dlmread('Team_82.dat');
A = 1 / (Cwall*R2);
B = 1 / (Cwall*R1);
C = 1 / (Cint*R1);
D = 1 / (Cint*Rwin);
% time dependent window 
tspan = linspace(0,1440,1441);  % calculate from t=0 up to t=1440
y0 = [18 20];    % initial conditions
% define your ode-function
odefunc = @odeFunction1;
[t,dTempdt] = ode45(odefunc, tspan, y0);

function dTempdt = odeFunction1(t,y)
% set of ordinary differential equations
% input:  A - constant
%         B - constant
%         C - constant
%         D - constant
%         t - the time variable
%         Tint - state variable
%         Twall - state variable
%output:  dTempdt - the set of equations
% initialize the set of equations
Tamb_t = interp1(tspan,Tamb,t);

Twall = y(1);
Tint = y(2);
dTempdt = zeros(2,1);
% define the set of equations
dTempdt(1) = A*Tamb_t - A*Twall + B.*Tint(:,1) - B.*Twall;
dTempdt(2) = C*Twall - C*Tint + D*Tamb_t - D*Tint;
end
end

[Hidde Kemperink]

Hi stefan!

Thanks for your answer but it is mandatory for us to use ode45 not ode1.

[Stephan]

it is usage of ode45 - ode1 is what i called my user defined function it. should i change the Name? no worries this is what you want

Answer 2

[Daniel]

Hi Hidde,

I composed one functioning script from all the comments in order to not confuse you.

Hope this works! ;)

Cheers,

Daniel

%% Import Data
opts = delimitedTextImportOptions("NumVariables", 1);
opts.DataLines = [1, Inf];
opts.Delimiter = ",";
opts.VariableNames = "Tamb";
opts.VariableTypes = "double";
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
Tamb = readtable("Team_82.dat", opts);
Tamb = table2array(Tamb);
clear opts

%% Solution

odefunc = @(t,y) odeFunction1(t,y); %output of odeFunction1 has to be dy/dt
tspan = linspace(0,length(Tamb),length(Tamb));
y0 = [18; 20];    % initial conditions

[t_solution,y_solution] = ode45(odefunc, tspan, y0);


function [dydt] = odeFunction1(t,y)
Tamb = evalin('base','Tamb');
R1 = 0.4;
R2 = 0.3;
Rwin = 0.5;
Cint = 1.0 * 10^4;
Cwall = 5.0 * 10^2;
A = 1 / (Cwall*R2);
B = 1 / (Cwall*R1);
C = 1 / (Cint*R1);
D = 1 / (Cint*Rwin);
Twall = y(1,1); %these are your state variables inside the state vector y
Tint = y(2,1); % you have two states (Twall and Tint)
tspan = linspace(0,length(Tamb),length(Tamb));
Tamb_t = interp1(tspan,Tamb,t);
dydt = zeros(2,1); %your output time derivative of the state vector
dydt(1,1) = A*Tamb_t - A*Twall + B*Tint - B*Twall;
dydt(2,1) = C*Twall - C*Tint + D*Tamb_t - D*Tint;
end

[Stephan]

Usage of global variables is not recomended: https://blogs.mathworks.com/videos/2010/03/08/top-10-matlab-code-practices-that-make-me-cry/

[Daniel]

You're right... dismiss the key word global then and use the following inside the ode-function:

Tamb = evalin('base','Tamb');

;)

[Stephan]

https://de.mathworks.com/help/optim/ug/passing-extra-parameters.html

[Daniel]

Stephan I agree with you. It's not the super cleanest way of coding but it gets the job done and isn't overloaded with complexity. Shouldn't confuse a first year student too much after all ;)

[Stephan]

I agree, but if we dont want to confse him, at first we should tell him that there is no inbuilt function ode_1

;-)

Answer 3

[PRAMOD KOTHMIRE]

Dear Hidde

I  want to plot t versus E where the formula for E is 
E= y(3)/int(y(3),t,0,2000)

Where to introduce the code for above formula in the following code so as to get the plot for t vs E.





function dydt = vdproc2(t,y)
global ntank Q V cainit tp
dydt=zeros(ntank,1)
cadum=y
if(t<tp)
cao=cainit
else
    cao=0
end
dydt(1)=Q/V*(cao-cadum(1))
for i=2:ntank
    dydt(i)=Q/V*(cadum(i-1)-cadum(i))
end
end

clear
clear all
clc

global ntank Q V cainit tp
tp=0.5;
Q=0.000006;
V=0.006;
tBegin = 0; % time begin
tEnd = 1000; % time end
Nio=0.00000452;
ntank=3;
cainit=Nio/(Q*tp)


[t,y] = ode45(@vdproc2,[tBegin tEnd],zeros(ntank,1));
 plot(t,y(:,3));
 
 ********************************************************************