MATLAB: Solving equations when variable is included in mu and sigma of distribution function

equationisalwayssymenginevariable in mu and sigmavpasolve

Hello
I'm trying to solve the following equation
c = 1.60
q = 0.08
equity = 0.08
p = 3.26
sigma = 0.1
mu = 0.08
    syms r
    vpasolve((c*q-equity)*(1+r)/p==(-p/(1+r)+p*(1-normcdf((c*q-equity)*(1+r)/p,(c*mu-equity)*(1+r)/p,c*(1+r)/p*sigma)))*(c*(1-normcdf(q,mu,sigma))-(c*q-equity)*normpdf(q,mu,sigma))/((p/(1+r))^2*normpdf(q,mu,sigma)-c^2*normpdf((c*q-equity)*(1+r)/p,(c*mu-equity)*(1+r)/p,c*(1+r)/p*sigma)),r,[0 1],'random',true)
However, Matlab returns an error message
symengine
Cannot prove '(1809252434946447*r)/36729926713742350 + 1809252434946447/36729926713742350 == 0 &
0.0000000000000000025671305227897154891565148192726*r + 0.0000000000000000025671305227897154891565148192726 <
- (12505220420059*r)/554413988131960 - 12505220420059/554413988131960' literally. To test the statement
mathematically, use isAlways.
   sym/subsindex (line 762)
                  X = find(mupadmex('symobj::logical',A.s,9)) - 1;
   sym/privsubsasgn (line 1031)
                  L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);
   sym/subsasgn (line 868)
              C = privsubsasgn(L,R,inds{:});
   normcdf>localnormcdf (line 100)
      p(sigma==0 & x<mu) = 0;
   normcdf (line 46)
  [varargout{1:max(1,nargout)}] = localnormcdf(uflag,x,varargin{:});
My guess is that Matlab is trying to test sigma==0 & x<mu, but it cannot because the symbolic variable r is included in sigma and mu. However, sigma is always >0 because of the range of r I have set.
How can I solve the equation in this case?

Best Answer

normpdf() and normcdf() are numeric routines. They are coded in a way that cannot succeed if the mu or sigma are symbolic. They internally need to test sigma == 0 & mu < 0 and also mu <= 0, which will fail for symbolic values that have insufficient assertions active (and even then, the assertions might not be enough.) It is okay to use symbolic values for the x for normpdf and normcdf though.
The above is why the vpasolve() is failing; it is failing while looking at the normpdf() expression before even getting into the vpasolve() call.
The solution is to turn it into a numeric system: transform vpasolve(A(r) == B(r),r) into @(r) A(r) - B(r) and apply a numeric solver to that.
Your vpasolve() is for r from 0 to 1. There is no solution to your system in that range. The solution is at approximately r = 1.008849984 .
Related Question