Numerical Solver in Matlab

The code below is part of a script that finds x to the nearest degree (BundleAngleArea), which lies between 0 and pi/2. The concept is as follows: A random 'BundleAngleArea' is generated. Integration of a certain expression from 0 to x is equal BundleAngleArea. The image below represents the area when x = pi/2 (BundleAngleArea = 0.5).

The expression contains three constants: BundleLength(25*10^(-3)), BundleWidth(8.4*10^(-3)) and FODParameter(0.5), which are decided and input by the user at the very beginning of the script.

As many 'BundleAngles' have to be generated (10000 - 100000 times), the expression was opened up to speed up Matlab solver and can be found in 'Constant' and 'd' of the code. However, I've encountered two problems:

i) I'll like to ask if there are any quicker methods of solving the equation other than the use of 'solve', as the speed remains very slow. I've attempted to use the function 'vpasolve' as well with no improvement in speed. I've also considered the use of a polynomial fit between BundleAngle and BundleAngleArea to replace the trigonometric expression and reduce complexity. However, different 'FODParameter's would require very different polynomial fits and thus require multiple 'for' loops to avoid poor accuracy. Hence, this alternative also seems unfeasible. Would there be any method to obtain a much quicker solution to this problem?

ii) The current solver returns a real and very small (10^(-11)) complex component (hence the use of 'real' function in the code). I've verified the results that come out of the real component and they seem to be fine. I'm quite curious as to why Matlab produces a complex part.

The code is as follows:

for BundleNumber = 1:10000

BundleAngleArea = rand/2;                                                              

Constant = 2*BundleAngleArea*((2/3)*(2*FODParameter - 1)*(BundleWidth - BundleLength) + (BundleLength + BundleWidth)) - BundleLength + (2*FODParameter - 1)*BundleLength*(2/3) ;

d = real(double(solve(BundleWidth*sin(x)- BundleLength*cos(x) + (2*FODParameter - 1)*(BundleLength*(cos(x) - (1/3)*cos(3*x)) + BundleWidth*((1/3)*sin(3*x) + sin(x))) - Constant == 0, x)));

BundleAngle(1,BundleNumber) = round(d(d>0)*180/pi);

end

Answer 1

[Roger Stafford]

As you have discovered the hard way, if you want efficient computation, it is not wise to use 'solve' for each of 10000 to 100000 different cases. That is an inherently slow process. Instead you should use 'solve' just once with general symbols for all four of your parameters, 'BundleAngleArea', 'BundleLength', 'BundleWidth', and 'FODParameter'. (I would suggest shorter symbols such as 'A', 'L', 'W', and 'F', respectively in order to simplify the solution from 'solve'.) On my ancient version of matlab the solution is of the form: twice the arctangent of the roots of a certain sixth-order polynomial whose seven coefficients are each simple expressions in terms of these four parameters. Therefore your code should simply evaluate these coefficients for the three fixed parameters and each of the varying BundleAngleArea's, and then call on matlab's numerical 'roots' function each time. The 'roots' function will produce six roots and your task will be to select the correct one of the six. Very likely only one of them will be real-valued and lying within the range 0 to 1 which would be necessary to produce an x value between 0 and pi/2. I would be greatly surprised if the time to compute seven coefficients, call on 'roots', and choose the proper root, isn't much, much less than a call on 'solve' with given numerical parameters.