Um, you don't need to do it that way? In fact, solving it that way is just a bit silly?
This problem reduces to something like
FCF(1)/(1+x) + FCF(2)/(1+x)^2 + ...
Essentially, this is a simple polynomial in 1/(1+x).
Just transform the problem using
Then your problem reduces to
FCF(1)*u + FCF(1)*u^2 + FCF(1)*u^3 + FCF(1)*u^4 + FCF(1)*u^5 + FCF(1)*u^6 - PV == 0
Solve it trivially as:
u = roots([flip(FCF),-PV]);
x = 1./u - 1
x =
-1.66107434285283 + 0i
-1.37767980674422 - 0.526238359452146i
-1.37767980674422 + 0.526238359452146i
-0.715412642705487 - 0.694273349944975i
-0.715412642705487 + 0.694273349944975i
0.0139259084189054 + 0i
Then pick out the root(s) that you care about. Here, it seems you want only the real root that lies in the interval [0,1] for x.
There is only one root in that interval.
x((imag(x) == 0) & (real(x) >= 0) & (real(x) <= 1))
ans =
0.0139259084189054
It happens to be the same root that vpasolve found, but in a tiny fraction of the time.
Now, while you cannot use roots in a vectorized form unless you use tools like cellfun, who cares? Just use a loop. Since a call to roots will be vastly more efficient than use of vpasolve, a loop is simple and efficient.
tic,u = roots([flip(FCF),-PV]);x = 1./u - 1;x=x((imag(x) == 0) & (real(x) >= 0) & (real(x) <= 1));toc
Elapsed time is 0.001027 seconds.
5000 such calls will take roughly 5 seconds, even on my old, slow CPU.
Compared to the same use of vpasolve as you wrote it, roots is roughly 175 times faster.
Best Answer