How to solve an equation system using a interp1-variable

7 views (last 30 days)
Anton Müller on 16 Sep 2021
Commented: Anton Müller on 17 Sep 2021
Hello everyone,
I'm trying to solve equations in a Matlab function that is part of a Simulink model.
One of those two equations is an interpolation of a given array.
% variable definiton
% T function argument
% p = f(T) values are given in form of a 69x2 double
% a | b | c | h constants / parameters (defined or calculated beforehand)
% equations
% eqn 1
% p = interp1(p(:,1),p(:,2),T)
% eqn 2
% h = a * T + p * (b + c * T)
% As I am used to solve equations with my pocket calculator I quickly found something like this while browsing:
syms T
solve (h == a * T + interp1(p(:,1),p(:,2),T) * (b + c * T),T);
% Unfortunately I do not have the required toolbox to solve equations which is why I need to find a work-around
% my next idea was to use an iteration to solve this (which I did not formulate and test until now)
% it would be something like that though:
for T = 1:100
p = interp1(p(:,1),p(:,2),T);
h_iter = a * T + p * (b + c * T);
if h_iter > h - h_marge | h_iter < h + h_marge
break
end
end
% with h_marge being a predefined value for accuracy tolerance
I would like to avoid such an iterational process since the iteration has to be done each simulation-step which could cause long simulation times
depending on the accuracy and the initial condition(s) I determine.
Thus, I wanted to ask you if there is another method / solution to solving that problem.
Look forward to hearing from you and thanks in advance.
Kind regards
Anton

John D'Errico on 16 Sep 2021
Edited: John D'Errico on 16 Sep 2021
interp1 is incompatible with symbolic parameters. You CANNOT use interp1 with solve. Period. So wherever you found that snippet of code, they were wrong.
You can use a tool like fzero with interp1. fzero and interp1 are compatible.
Next, looking at your problem, this is actually a TWO variable problem, with two equations. That is, we have known constants, a,b,c,h. And there are two unknowns. Sadly, you seem to be using p in two places at once, both as an array of known elements, and as a variable. And that is confusing as hell.
So I'll make the assumption that you have two variales, I'll call then T and p_T. That is p_T is given as:
p_T = interp1(p(:,1),p(:,2),T)
So, if I knew the value of T, then I could compute P_T. And then we have a second equation:
h = a * T + p_T * (b + c * T)
Now we can just substitute P_T in there. This reduces the problem to ONE equation, in one unknown, where you can use fzero. I've re-arranged it to be something you can see as a problem for fzero.
a * T + interp1(p(:,1),p(:,2),T) * (b + c * T) - h = 0
Now, given the variables a,b,c,h and the array p are defined as MATLAB variables in your workspace, set it up as:
fun = @(T) a * T + interp1(p(:,1),p(:,2),T) * (b + c * T) - h;
Tfinal = fzero(fun,[p(1,1),p(end,1)]);
That call will set a bracket on T, such that T must always lie between P(1,1) and P(69,1). That must make sense, since interp1 cannot intelligently evaluate its function beyond those limits anyway.
Again, a,b,c,h and p must all be previously defined in your workspace for that to work.
Anton Müller on 17 Sep 2021
Hello John,
thank you very much for your fast response. I simplified the (originally thermodynamical) problem into the form of my question. While doing that I overlooked that I used p twice. You're completely correct they are two seperate values:
p = [69x2] double
p_T
Your proposed solution worked out fine in an example but after transforming it back to the original problem I got the following error:
Error using fzero (line 274)
Function values at the interval endpoints must differ in sign.
Error in Test (line 80)
Tfinal = fzero(fun,[p(1,1),p(end,1)]);
I found the solution after reading a bit about the fzero command, though. If it is of any interest:
eqn2: h = a * T + p_T * (b + c * T)
The p_T part would correctly be:
d * (p_T /(p_b - p_T))
d and p_b are constants defined prior. This would mean the correct expression eqn2' for eqn2 would be
eqn2': h = a * T + d * (p_T /(p_b - p_T)) * (b + c * T)
The problem was that p(end,1) > p_b. Thus, it could not determine whether there are zeros in the interval. I could adjust to that by deleting p(69,1) as it wasn't relevant for the application.
Thank you for your time and effort.
Kind regards
Anton

R2021a

Community Treasure Hunt

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

Start Hunting!