Using Ode45 to solve differential equation with time dependent variable

107 Ansichten (letzte 30 Tage)
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!
  8 Kommentare
Hidde Kemperink
Hidde Kemperink am 25 Okt. 2019
Bearbeitet: Hidde Kemperink am 25 Okt. 2019
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
Stephan am 25 Okt. 2019
Bearbeitet: Stephan am 25 Okt. 2019
I changed the name - now it should be clear that i use ode45...

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Stephan
Stephan am 25 Okt. 2019
Bearbeitet: Stephan am 25 Okt. 2019
% 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
  2 Kommentare
Hidde Kemperink
Hidde Kemperink am 25 Okt. 2019
Hi stefan!
Thanks for your answer but it is mandatory for us to use ode45 not ode1.
Stephan
Stephan am 25 Okt. 2019
Bearbeitet: Stephan am 25 Okt. 2019
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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (2)

Daniel
Daniel am 25 Okt. 2019
Bearbeitet: Daniel am 25 Okt. 2019
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
  5 Kommentare
Daniel
Daniel am 25 Okt. 2019
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
Stephan am 25 Okt. 2019
I agree, but if we dont want to confse him, at first we should tell him that there is no inbuilt function ode_1
;-)

Melden Sie sich an, um zu kommentieren.


PRAMOD KOTHMIRE
PRAMOD KOTHMIRE am 12 Mai 2020
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));
********************************************************************

Kategorien

Mehr zu Startup and Shutdown finden Sie in Help Center und File Exchange

Produkte


Version

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by