4 views (last 30 days)

Show older comments

I wrote a problem set for first time matlab programmers, which involved using the solve routine. Unfortunately, they discovered the following problem, which did not endear them to either to me or matlab.

In R2016a, the solve command has the newly acquired facility to generate a numeric solution when it can't find a symbolic one, which is a very nice feature, except that it appears to be extremely fragile, treating essentially identical problems quite differently. In the example below, we compute a derivative manually, and also using matlab's symbolic diff routine. The two derivatives are identical except for rounding error of order 1e-16. In one case, the solve routine computes the correct numeric solution. In the other, it generates 17 symbolic solutions, which appear to us to be nonsense. Is this intentional behavior? Or a bug?

%%set parameters

xunder = 3;

T = 15;

b = 0.9;

r = 0.2;

rho = 2.125;

kappa = 2;

syms x nu t;

%%create base functions

c_t = @(nu,t) nu * ((1+r) * b)^t * (1+r) * (1-b);

x_t = @(nu,t) (1+r)^t * (xunder - nu*(1-b^t));

price = @(x) (x/kappa)^(-1/rho);

scrapVal = @(x) price(x) * x;

%%differentiate the scrapValue function

% using matlab

scrapPrimeMatlab = matlabFunction(diff(scrapVal(x)),'Vars',x);

% by hand

scrapPrimeManual = @(x) (1 - 1/rho)*(kappa./x).^(1/rho);

% show they are essentially equal

scrapPrimeDiff = @(x) scrapPrimeMatlab(x) - scrapPrimeManual(x);

X = linspace(2,10,100);

disp(['Maximum difference between matlab''s deriv and manual deriv is ' num2str(max(abs(scrapPrimeDiff(X)))) ]);

disp(' ');

%%Create Control/State Functions

eqControlManual = @(x,nu) b * scrapPrimeManual(x) - 1/(c_t(nu,T));

eqControlMatlab = @(x,nu) b * scrapPrimeMatlab(x) - 1/(c_t(nu,T));

eqState = @(x,nu) x_t(nu,T+1)-x;

% output solutions

solnManual = solve([eqControlManual(x,nu)==0, eqState(x,nu)==0]);

disp(['Using the manual derivative, the solve command generates the numeric solution: x=' num2str(eval(solnManual.x)) '; nu=' num2str(eval(solnManual.nu)) ]);

disp(' ');

disp('Now running the solve command, replacing the manual derivative with matlab''s derivative')

solnMatlab = solve([eqControlMatlab(x,nu)==0, eqState(x,nu)==0]);

solnMatlab

% different solutions, same objective function only differ by taking

% derivative by hand with scrapPrimeManual and with matlabFunction(diff(.))

% where the derivative by hand and matlab differentiation are essentially

% equal

John D'Errico
on 3 Jun 2016

Edited: John D'Errico
on 3 Jun 2016

Not a bug. Part of your problem is in misunderstanding floating point arithmetic, where the number 0.9 is NOT represented exactly. Nor, for example, is the number 2.5e1./1.7e1.

So lets look to see where the difference arises.

scrapPrimeMatlab

scrapPrimeMatlab =

@(x)x.*1.0./(x.*(1.0./2.0)).^(2.5e1./1.7e1).*(-4.0./1.7e1)+1.0./(x.*(1.0./2.0)).^(8.0./1.7e1)

scrapPrimeManual

scrapPrimeManual =

@(x)(1-1/rho)*(kappa./x).^(1/rho)

These are very different expressions. Subtle differences in those exponents down in the 16th decimal digit will be huge in terms of being able to solve the equations.

The fact is, symbolic algebra is not always intelligent algebra. Very often there will be terms that should combine or subtract out, but subtle differences down in those least significant bits will cause problems.

You would probably have been better off by defining numbers like 0.9 as a symbolic 0.9. Then MATLAB will be smart enough to know the number is EXACTLY 9/10. Too many of your computations were done in floating point arithmetic, and only at the very end do you tell MATLAB that this is intended to be a symbolic problem.

For example, suppose I tried this:

sym(0.9 - 0.3 - 0.6)

ans =

1/9007199254740992

Should you be surprised? Not really, IF you understand how a tool like MATLAB must do that computation. First it computes the result

0.9 - 0.3 - 0.6

ans =

1.1102e-16

as a FLOATING point number. Then it passes that into sym. Oh, by the way, I wanted you to make that result symbolic. MATLAB cannot read your mind, at least not yet.

Sorry, but virtually computer code can be screwed up by someone who is sufficiently diligent at the task.

John D'Errico
on 3 Jun 2016

Hi Derek,

It is good that you are asking questions, wanting to understand.

I think your first problem MAY be in the use of functions too much here. Expressions that employ symbolic variables are already functions. For example, what if we just use syms all the way? Integers we can leave alone.

xunder = 3;

T = 15;

b = sym('0.9');

r = sym('0.2');

Note that 2.125 is indeed exactly representable as a double, so I could probably get away with leaving r as a double. But it is better to be safe than sorry.

rho = sym('2.125');

kappa = 2;

syms x nu t;

%%create base functions

c_t = nu * ((1+r) * b)^t * (1+r) * (1-b);

x_t = (1+r)^t * (xunder - nu*(1-b^t));

price = (x/kappa)^(-1/rho);

scrapVal = price * x;

Lets see where that got us:

scrapVal

scrapVal =

x/(x/2)^0.47058823529411764705882352941176

Now, lets differentiate. (Computers can be dumb, so don't be surprised at this next pair of results.)

scrapPrime = diff(simplify(scrapVal))

scrapPrime =

0.73359229710742684336247382121409/x^0.47058823529411764705882352941176

scrapPrime2 = simplify(diff(scrapVal))

scrapPrime2 =

1.3856743389806951485735616622933/x^0.47058823529411764705882352941176 - 0.65208204187326830521108784107919/x^0.47058823529411764705882352941176

Yes, you and I can see the two results are the same. And with some effort, I can probably get MATLAB to combine the latter form. Not worth the brain sweat though. :)

eqControl = b*scrapPrime - 1/subs(c_t,t,T)

eqControl =

0.66023306739668415902622643909268/x^0.47058823529411764705882352941176 - 2.6270142080490851321549115216431/nu

eqState = subs(x_t,t,T+1)-x

eqState =

55.4652776685109248 - x - 15.062483246169511161754270763581*nu

Now just use solve.

soln = solve([eqControl==0, eqState==0])

soln =

nu: [1x1 sym]

x: [1x1 sym]

soln.nu

ans =

3.6277902362881364768614086810024

soln.x

ans =

0.82174801380353695115956262208444

In fact, I'd expect there to be no analytical solution. Pseudo-polynomial problems with fractional exponents like this almost never have an analytical solution.

Long ago I learned to go slow. Don't just throw a large block of code together and expect it to work and to do exactly as you wish, without thought on your part. Check EVERY line. Here, I don't just mean that it is mathematically correct. Does the result make sense? If symbolic,is it as simple as possible? (At least do so until you have the experience to know what you can get away with.) The point is, check each result for sanity.(I forgot to do that far too often in my early years. Ok, one truly notable time was enough to teach me that lesson.) Look to see if you could do something to improve it. For example, remember the subtle difference between simplify(diff)) and diff(simplify()). Computers can be just plain dumb sometimes. Yeah, some of the time, it is me being dumb. But not always.

As I wrote it above, see that I never had to resort to function handles. Those expressions are already functions of symbolic variables. But as I wrote them, see that by living as much as possible in the symbolic domain, MATLAB was happier. A happy computer is a good computer.

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

Start Hunting!