How to solve an equation system using a interp1-variable

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

Answer 1

[John D'Errico]

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]

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