Using interp1 within vpasolve

I'm having difficulty with the following situation trying to numerically solve an equation (with vpasolve, or something similar) that uses interp1-interpolated data within the solve. Consider:

with implementation:

syms Z;
vpasolve(Z.^4 == ( muL .* interp1(X,GH,Z) - muH .* interp1(X,GL,Z) )./( muL .* interp1(X,FH,Z) ...
    - muH .* interp1(X,FL,Z) ), Z);

where and are constants, so that the μ-values are constant, and F and G are 2D matrices squeezed to 1D. These 1D matrices are known on a grid Z = 5:15. Then, I'm attempting to solve this equation for Z, which presumably has real-numbered values as solution. I assumed vpasolve worked iteratively, with an initial guessed query of Z, such that interp1 with (syms Z) is valid to use inside vpasolve. Of course, this doesn't work.

Is there any way to work this? I've simplified things for ease of explanation, but suffice it to say that just doing a polynomial fit on F and G instead of using interp1 is not on the table: the data in the F and G grids need to only be interpolated for non-integer values of Z. Here's some data, just to make it easier to run if needed.

muH = 0.178125;
muL = 0.142500;
GL = [0.2684 0.2690 0.2696 0.2704 0.2714 0.2724 0.2734 0.2746 0.2756 0.2768 0.2781];
GH = [0.3258 0.3280 0.3307 0.3342 0.3381 0.3431 0.3474 0.3520 0.3568 0.3618 0.3670];
FL = 1E-6 .* [0.4234 0.4180 0.4100 0.4017 0.3932 0.3840 0.3766 0.3698 0.3628 0.3554 0.3480];
FH = 1E-6 .* [0.4234 0.4180 0.4100 0.4017 0.3932 0.3840 0.3766 0.3698 0.3628 0.3554 0.3480];

Answer 1

[Walter Roberson]

Turn it into a numeric system

z0 = 1;

muH = 0.178125;
muL = 0.142500;
GL = [0.2684 0.2690 0.2696 0.2704 0.2714 0.2724 0.2734 0.2746 0.2756 0.2768 0.2781];
GH = [0.3258 0.3280 0.3307 0.3342 0.3381 0.3431 0.3474 0.3520 0.3568 0.3618 0.3670];
FL = 1E-6 .* [0.4234 0.4180 0.4100 0.4017 0.3932 0.3840 0.3766 0.3698 0.3628 0.3554 0.3480];
FH = 1E-6 .* [0.4234 0.4180 0.4100 0.4017 0.3932 0.3840 0.3766 0.3698 0.3628 0.3554 0.3480];

best_z = fsolve(@(Z) Z.^4 - (( muL .* interp1(X,GH,Z,'linear','extrap') - muH .* interp1(X,GL,Z,'linear','extrap') )./( muL .* interp1(X,FH,Z,'linear','extrap') ...
    - muH .* interp1(X,FL,Z,'linear','extrap') )), z0);

Unrecognized function or variable 'X'.

Error in solution (line 10)
best_z = fsolve(@(Z) Z.^4 - (( muL .* interp1(X,GH,Z) - muH .* interp1(X,GL,Z) )./( muL .* interp1(X,FH,Z) ...

Error in fsolve (line 260)
fuser = feval(funfcn{3},x,varargin{:});

Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.

[Walter Roberson]

vpasolve() does not just proceed iteratively. It first checks to see whether the system is polynomial, and if so then uses solve() to find the roots and then uses vpa() on them. If not then it does a modified Newton-Raphson -- which requires that it be able to take the derivative of the expression.

Either way, remember that any parameter you give to a function is going to be fully evaluated before it has passed into the function. If you had, for example, sin(0.183) in the middle of the expression, then that would be evaluated to a double precision number and incorporated into the larger expression. What you pass as a symbolic expression to vpasolve() does not have "delayed evaluation".

It happens that for symbolic variable Z that Z.^4 evaluates out to something symbolic, and that symbolic result becomes part of the value that will get passed into vpasolve() -- but not until you get full evaluation of the expression. Z^4 - Z^5/Z + Z will get simplified away to Z, not taken intact into vpasolve()

This matters because you have interp1(X,GH,Z) in there, but interp1() is not defined for symbolic values -- and at the first pass when the expression is being evaluated before being passed in to vpasolve(), that matters. You cannot get around it by coding interp1(X,GH,double(Z)) because the Z would be an unresolved symbolic variable at that point, and double() would error on unresolved symbolic variables.

You would have to change your coding from symbolic expression to function handle, in order to delay evaluation of Z until some actual numeric value was passed in.

You can pass a function handle to vpasolve(). However, vpasolve() will try to transform the function handle into a symbolic expression by passing a symbolic variable into the function, and then you would have the same problem all over again with interp1()

... So you need to give up on the idea of using vpasolve(), unless you can code a symbolic alternative to interp1.

It is possible to build up a symbolic alternative to interp1() using piecewise(). Whether that is a good idea is a different matter.

[Adam Zieser]

Unfortunate :( Is there any way to numerically solve an equation using interpolated data from a grid in the desired way? Perhaps without vpasolve? Mathematica has no issue with this task in the slightest, using FindRoot[ ] with user-defined functions based on Interpolation[ ] of F & G functions.

[Walter Roberson]

The code I posted above will do a numeric interpolation, provided that you first define X.

[Adam Zieser]

Ah, yes, had to download the Optimization Toolbox. Thanks a bunch!

[Walter Roberson]

Could also have used fzero in this case.