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.
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];
采纳的回答
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);
5 个评论
Walter Roberson
2021-7-14
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
2021-7-14
Walter Roberson
2021-7-14
The code I posted above will do a numeric interpolation, provided that you first define X.
Adam Zieser
2021-7-14
Walter Roberson
2021-7-14
Could also have used fzero in this case.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Equation Solving 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
