matlab ode45 wont accept new parameter

% Set extra parameters
tf = 72; % h
h0 = 2; % m
g = 9.81;

%prompt = 'What is the leakage value? ';
%a = input(prompt);
%L = a*g*h;

% Solve the DAE
opts = odeset('Reltol',1e-5);
[t,h] = ode45(@tankf1,[0 tf],h0,opts);
% Plot the graph
plot(t,h)
title('Tank model solved by Method 1')
xlabel('Time (h)')
ylabel('Level (m)')
grid on

function P1 = P1function(t)
% P1FUNCTION: Upstream pressure (kPa) as a function of time (h)
if t<=10
    P1 = 139.5;
elseif t<=11
    P1 = 139.5 + 60.5*(t-10);
elseif t<=40
    P1 = 200;
elseif t<=55
    P1 = 200 - 30*(t-40)/15;
else  % t>55
    P1 = 170;
end
end

function dhdt = tankf1(t,h)
% TANKF1: Tank level model - Derivative function for solution method 1
% (substitution of algebraic variables)
% Tank model parameters and constants
A = 12.5; % m^2
CV1 = 3.41; % (m^3/h)/kPa^0.5
CV2 = 3.41; % (m^3/h)/kPa^0.5
P0 = 100; % kPa
P3 = 100; % kPa
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
P1 = P1function(t); % kPa

% Algebraic equations for F1, F2 and P2 have been substituted into the
% differential equation for dh/dt to create a pure ODE system:
dhdt = 1/A*(CV1*sqrt(P1-P0-rho*g*h/1000) - CV2*sqrt(P0+rho*g*h/1000-P3))-L;
end

I have the code above shown, and it is a simulation of a tank with water levels at a certain time, up to 72 hours. The issue im having is I need to simualte it with a leakage value, this is calculated as shown at the top of code using a*g*h. The issue im having it when I try to put it in the dhdt equation in the function for tankf1, it keeps throwing an ode error and cant figure out how to alter this. NB I have to use the ode to solve the equation and draw the graph.

Answer 1

[madhan ravi]

EDITED

If you want to pass a new parameter to a function parametrize your function

% Set extra parameters
tf = 72; % h
h0 = 2; % m
g = 9.81;
a=input('leakage value : ');
L=a*g*h;
%prompt = 'What is the leakage value? ';
%a = input(prompt);
%L = a*g*h;
% Solve the DAE
opts = odeset('Reltol',1e-5);
[t,h] = ode45(@(t,h)tankf1(t,h,L),[0 tf],h0,opts);
% Plot the graph
plot(t,h)
title('Tank model solved by Method 1')
xlabel('Time (h)')
ylabel('Level (m)')
grid on
function dhdt = tankf1(t,h,L)
% TANKF1: Tank level model - Derivative function for solution method 1
% (substitution of algebraic variables)
% Tank model parameters and constants
A = 12.5; % m^2
CV1 = 3.41; % (m^3/h)/kPa^0.5
CV2 = 3.41; % (m^3/h)/kPa^0.5
P0 = 100; % kPa
P3 = 100; % kPa
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
P1 = P1function(t); % kPa
% Algebraic equations for F1, F2 and P2 have been substituted into the
% differential equation for dh/dt to create a pure ODE system:
dhdt = 1/A*(CV1*sqrt(P1-P0-rho*g*h/1000) - CV2*sqrt(P0+rho*g*h/1000-P3))-L;
end

function P1 = P1function(t)         % should be defined at last
% P1FUNCTION: Upstream pressure (kPa) as a function of time (h)
if t<=10
    P1 = 139.5;
elseif t<=11
    P1 = 139.5 + 60.5*(t-10);
elseif t<=40
    P1 = 200;
elseif t<=55
    P1 = 200 - 30*(t-40)/15;
else  % t>55
    P1 = 170;
end
end

example of the graph produced:

[Mitul Dattani]

The issue I had with putting it in the function was I need a to be an inputted value but when i add the input method into the function it just keeps repeating and never ends

[Mitul Dattani]

function dhdt = tankf1(t,h)
% TANKF1: Tank level model - Derivative function for solution method 1
% (substitution of algebraic variables)
% Tank model parameters and constants
A = 12.5; % m^2
CV1 = 3.41; % (m^3/h)/kPa^0.5
CV2 = 3.41; % (m^3/h)/kPa^0.5
P0 = 100; % kPa
P3 = 100; % kPa
rho = 1000; % kg/m^3
g = 9.81; % m/s^2

input = 'What is the leakage value? ';
a = input(input);
L=a*g*h;         % change L to your value

P1 = P1function(t); % kPa

% Algebraic equations for F1, F2 and P2 have been substituted into the
% differential equation for dh/dt to create a pure ODE system:
dhdt = 1/A*(CV1*sqrt(P1-P0-rho*g*h/1000) - CV2*sqrt(P0+rho*g*h/1000-P3))-L;
end

This is initally what I had tried, but it gives me the error:

Index exceeds matrix dimensions.

Error in Method1>tankf1 (line 31)
a = input(input);

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 Method1 (line 10)
[t,h] = ode45(@tankf1,[0 tf],h0,opts);

[madhan ravi]

wierd , it worked fine for me , ok save the corresponding functions with corresponding name and run the main code , or paste the full code what iu are trying here

[Mitul Dattani]

% Set extra parameters
tf = 72; % h
h0 = 2; % m
g = 9.81;

% Solve the DAE
opts = odeset('Reltol',1e-5);
[t,h] = ode45(@tankf1,[0 tf],h0,opts);
% Plot the graph
plot(t,h)
title('Tank model solved by Method 1')
xlabel('Time (h)')
ylabel('Level (m)')
grid on
function P1 = P1function(t)
% P1FUNCTION: Upstream pressure (kPa) as a function of time (h)
if t<=10
    P1 = 139.5;
elseif t<=11
    P1 = 139.5 + 60.5*(t-10);
elseif t<=40
    P1 = 200;
elseif t<=55
    P1 = 200 - 30*(t-40)/15;
else  % t>55
    P1 = 170;
end
end
function dhdt = tankf1(t,h)
% TANKF1: Tank level model - Derivative function for solution method 1
% (substitution of algebraic variables)
% Tank model parameters and constants
A = 12.5; % m^2
CV1 = 3.41; % (m^3/h)/kPa^0.5
CV2 = 3.41; % (m^3/h)/kPa^0.5
P0 = 100; % kPa
P3 = 100; % kPa
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
P1 = P1function(t); % kPa
% Algebraic equations for F1, F2 and P2 have been substituted into the
% differential equation for dh/dt to create a pure ODE system:
dhdt = 1/A*(CV1*sqrt(P1-P0-rho*g*h/1000) - CV2*sqrt(P0+rho*g*h/1000-P3));
end

So this is the full code which generates a graph, which is a sim of a tank and its water level based on pressure (in simplistic terms). There is a leakage value between 0.01 and 0.2 which is user inputted to try different values up to 0.2. So if i add it into the function as a user inputted value it throws me that error. where L = a*g*h, and the new dhdt should be equal to dhdt - L.

[madhan ravi]

just see my edited answer , copy it , paste it in a script , save it and just run it

[Mitul Dattani]

Error using odearguments (line 95)
@(T,H)TANKF1(T,H,L) returns a vector of length 529, but the length of initial conditions vector is 1. The
vector returned by @(T,H)TANKF1(T,H,L) and the initial conditions vector must have the same number of
elements.

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

Error in Method1 (line 12)
[t,h] = ode45(@(t,h)tankf1(t,h,L),[0 tf],h0,opts);

gets me this error already tried that

[Mitul Dattani]

The problem im having is wth the ode the input for L i can arrange later its just the reason as to why this error message is appearing

[madhan ravi]

Yes even I am having difficulties in figuring it out because h is returned by ode solver

[Mitul Dattani]

yea thats the issue but cant see a way to use L to change the graph

[Steven Lord]

Mitul,

When you wrote this:

input = 'What is the leakage value? ';
a = input(input);
L=a*g*h;         % change L to your value

there is one problem with the code and one potential problem with the workflow.

Problem with the code: because you've defined a variable named input, you cannot call the input function. Rename the variable input to something like promptString.

Potential problem with the workflow: as written each and every time ode45 calls your ODE function tankf1 it will prompt the user to enter a value. Since the ODE solver will call your ODE function many times, this seems inefficient. Unless the leakage value changes over time in a way that's not computable from the previous leakage value and the current time, I'd prompt the user once for the value before you call ode45 and pass that value into your ODE function as a parameter as madhan demonstrated.

[madhan ravi]

Thank you very much Steven Lord didn't notice the variable input ;-)