Does MATLAB have a feature to show step-by-step solution?

Hello everyone.

I have a pretty complicated equation where I need to find the solution for B. I substituted the values of variables and separated the equation into 4 big "parts" to make it easier to read.

100-B==
100*exp(-(0.05)*m*(0.15))
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)-(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))
-exp(-((0.03)-(0.01))*m*(0.15))*B
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)+(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))
+(0.15)*(0.05)*100*symsum(exp(-(0.05)*(m-j)*(0.15)) ...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)-(0.2)^2/2)*(m-j)*(0.15)) ...
    /((0.2)*sqrt((m-j)*(0.15))))),j,0,m-1)
-(0.15)*((0.03)-(0.01))*exp((0.01)*m*(0.15))*B*symsum(exp(-((0.03)+(0.01))*(m-j)*(0.15))* ...
    normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)+(0.2)^2/2)*(m-j)*(0.15)) ...
    /((0.2)*sqrt((m-j)*(0.15))))),j,0,m-1)

where m denotes time and B denotes "price". I used vpasolve(<equation>,B) to solve the equation.

When I used m=4, it returns a real value for the solution, but starting from m=5, an imaginary component showed up. B is not supposed to have imaginary values at any point of time m.

I tried decomposing the equation and evaluating the value at each "part" at m=5 to figure out where the imaginary component originated from, but I cannot seem to find anything that results in a negative value inside the square root. Does MATLAB have a feature to show step-by-step solution?

Thank you very much!

[Torsten]

If vpasolve iterates B such that it becomes negative, log(B) will become complex-valued.

Restrict B to be positive:

assume(B,'positive') 

[Janice]

Thank you for the quick reply!

I tried using the positive assumption, but unfortunately it still has the same output..

[Christine Tobler]

When you specify the symbolic variable B, you can use

syms B real

instead of just

syms B

That should prevent returning complex values - although there's a chance that it would just return no solution if it otherwise finds only complex solutions.

By the way, what is j here? Is it another constant defined before this line is run, or does it denote sqrt(-1)?

[Janice]

Thank you very much for the reply!

I have tried using

syms B real

But somehow it still returns complex solution..

Sorry I forgot to mention that, but j is just a symbol for the symsum function where it ranges from to (m-1). Could there be a possibility that j affects the existence of a complex solution, even though it is already defined as a symbol?

[Christine Tobler]

If it's already defined as a symbol, j won't have any effect like that. I was just wondering as an outside possibility - if you had forgotten to define j, it might have been that its default value of sqrt(-1) was used instead.

[Janice]

Oh I see, thank you very much!

Answer 1

[Torsten]

m = 5;
B0 = 1;
B = fsolve(@(B)fun(B,m),B0)

B = 1

B = 1.0000

B = 2

B = 2.0000

B = 4.5000

B = 4.5000

B = 10.7500

B = 10.7500

B = 26.3750

B = 26.3750

B = 65.4375

B = 65.4375

B = 126.9775

B = 80.8225

B = 80.8225

B = 81.8887

B = 81.8887

B = 81.8346

B = 81.8346

B = 81.8344

B = 81.8344

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

B = 81.8344

fun(B,m)

B = 81.8344

ans = 1.9267e-10

function res = fun(B,m)
B
sum1 = 0.0;
sum2 = 0.0;
for j=0:m-1
    sum1 = sum1 + exp(-(0.05)*(m-j)*(0.15)) ...
        *normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)-(0.2)^2/2)*(m-j)*(0.15)) ...
        /((0.2)*sqrt((m-j)*(0.15)))));
    sum2 = sum2 + exp(-((0.03)+(0.01))*(m-j)*(0.15))* ...
        normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)+(0.2)^2/2)*(m-j)*(0.15)) ...
        /((0.2)*sqrt((m-j)*(0.15)))));
end
res = 100*exp(-(0.05)*m*(0.15))...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)-(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))...
    -exp(-((0.03)-(0.01))*m*(0.15))*B...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)+(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))...
    +(0.15)*(0.05)*100*sum1 ...
    -(0.15)*((0.03)-(0.01))*exp((0.01)*m*(0.15))*B*sum2 - (100-B);
end

[Janice]

Thank you very much!

I'm really sorry I am fairly new to coding, but I am having trouble with defining the function.

I tried running the code

function res = fun(B,m)
B
sum1 = 0.0;
sum2 = 0.0;
for j=0:m-1
    sum1 = sum1 + exp(-(0.05)*(m-j)*(0.15)) ...
        *normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)-(0.2)^2/2)*(m-j)*(0.15)) ...
        /((0.2)*sqrt((m-j)*(0.15)))));
    sum2 = sum2 + exp(-((0.03)+(0.01))*(m-j)*(0.15))* ...
        normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)+(0.2)^2/2)*(m-j)*(0.15)) ...
        /((0.2)*sqrt((m-j)*(0.15)))));
end
res = 100*exp(-(0.05)*m*(0.15))...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)-(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))...
    -exp(-((0.03)-(0.01))*m*(0.15))*B...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)+(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))...
    +(0.15)*(0.05)*100*sum1 ...
    -(0.15)*((0.03)-(0.01))*exp((0.01)*m*(0.15))*B*sum2 - (100-B);
end

but I get an error message saying

"Error: Function definition are not supported in this context. Functions can only be created as local or nested functions in code files."

I am not quite sure where I did wrong..

[Image Analyst]

Was that function in its own m-file called fun.m? Or was it part of another script or m-file? Sounds like the latter. If you have a script and functions in the same m-file, the functions must all end with an "end" statement and come AFTER the script.

[Torsten]

Load all this in the editor and run it.

By "run" I mean click on the green RUN arrow to the right at the task bar.

m = 1:100;
B0 = 1;
B = arrayfun(@(m)fsolve(@(B)fun(B,m),B0),m)
%fun(B,m)
plot(m,B)

function res = fun(B,m)
sum1 = 0.0;
sum2 = 0.0;
for j=0:m-1
    sum1 = sum1 + exp(-(0.05)*(m-j)*(0.15)) ...
        *normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)-(0.2)^2/2)*(m-j)*(0.15)) ...
        /((0.2)*sqrt((m-j)*(0.15)))));
    sum2 = sum2 + exp(-((0.03)+(0.01))*(m-j)*(0.15))* ...
        normcdf(-((log(B)+(0.01)*m*(0.15)-(log(B)+(0.01)*(m-j)*(0.15))+((0.05)-(0.03)+(0.2)^2/2)*(m-j)*(0.15)) ...
        /((0.2)*sqrt((m-j)*(0.15)))));
end
res = 100*exp(-(0.05)*m*(0.15))...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)-(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))...
    -exp(-((0.03)-(0.01))*m*(0.15))*B...
    *normcdf(-((log(B)+(0.01)*m*(0.15)-log(100)+((0.05)-(0.03)+(0.2)^2/2)*m*(0.15))/((0.2)*sqrt(m*(0.15)))))...
    +(0.15)*(0.05)*100*sum1 ...
    -(0.15)*((0.03)-(0.01))*exp((0.01)*m*(0.15))*B*sum2 - (100-B);
end

[Janice]

I just ran the code and it worked perfectly.

Thank you very much for all your help!