# MATLAB giving me a 16x1 complex double answer?

2 views (last 30 days)
Alma Morgan on 24 Feb 2023
Answered: Walter Roberson on 31 May 2023
I have two non-equations and two unknowns. I could like to be able to solve for the varibles. They sould be a simple answers, like 100 and 75, for examples.
I'm getting 16x1 double complex for an answer, and I'm not even sure what to make of it.
syms w h
L = 101.6 %mm
T = 8.7867e+06 %N*mm
P = 89.81 %N
G = 26000 %MPa / N/mm
E = 68900 %MPa / N/mm
eq1 = @(w,h) (T*L)/[w*h^3*[(16/3)-3.36*(1-(h^4/(12*w^4)))]*G] == 0.1
eq2 = @(w,h) (P*L^3)/(48*E*(w*h^3)/12) == 1
[val,sol] = solve({eq1,eq2},[w,h])
vals = double(val)
vals2 = double(sol)

Torsten on 24 Feb 2023
Replace
syms w h
by
syms w h real positive

Walter Roberson on 31 May 2023
You are using floating point inputs, each of which is inexact. For example is 89.81 intended to represent exactly 8981/100 or is it 1579954228649001/17592186044416 (as MATLAB would create it to be), or is it intended to represent "some number not precisely known that is between 89+805/1000 inclusive and 89+815/1000 exclusive (which is the interpretation if it is scientific notation.)
You then use solve(), but solve() is for the case where you are trying to calculate indefinitely precise results in closed-form formulas -- for example solve would generate (sqrt(5)-1)/2 instead of 0.618033988749895 if appropriate .
You then use double() -- suggesting that you do not care about the full precision, and only wanted an approximation after all.
If you have numbers that are not exact to start with, it is a mistmatch to try to find indefinitely correct solutions... especially if you discard the accuracy afterwards.
What to do then? Use vpasolve(). And throw away the results you do not want
syms w h
L = 101.6 %mm
L = 101.6000
T = 8.7867e+06 %N*mm
T = 8786700
P = 89.81 %N
P = 89.8100
G = 26000 %MPa / N/mm
G = 26000
E = 68900 %MPa / N/mm
E = 68900
eq1 = @(w,h) (T*L)/[w*h^3*[(16/3)-3.36*(1-(h^4/(12*w^4)))]*G] == 0.1
eq1 = function_handle with value:
@(w,h)(T*L)/[w*h^3*[(16/3)-3.36*(1-(h^4/(12*w^4)))]*G]==0.1
eq2 = @(w,h) (P*L^3)/(48*E*(w*h^3)/12) == 1
eq2 = function_handle with value:
@(w,h)(P*L^3)/(48*E*(w*h^3)/12)==1
[val,sol] = vpasolve([eq1(w,h),eq2(w,h)],[w,h]);
mask = imag(val) == 0 & imag(sol) == 0;
vals = 2×1
-0.9269 0.9269
vals2 = 2×1
-7.1706 7.1706