Matlab Solve and simplify functions not working

I'm using the 'solve' function to solve a system of equations. It is solving them but not simplifying. That's when I have tried to use the 'simplify' function. Simplifying is not working. It stays the same as before instead of being solved. If I copy/paste the expression, it solves it though. Any suggestions for getting it to simplify to one number without having to manually copy/paste for each one? For example, f21_sol(1) = -(760*(348475767795489804396578940311480912*3^(1/2) + 168032379110159424438402050627659277))/(3*(3858879662023561551604431872957445612*3^(1/2) - 3892977404894321714789280647602584263)), which is absolutely solvable for a single number.

Answer 1

[Paul]

Hi Brittany,

Run the code:

foot_forces
f21_sol

f21_sol = 

One option is to use vpa to convert to a VPA number. However, the result is still a sym, not a numeric that you could use in other Matlab functions. See the linked doc page for more information.

vpa(f21_sol,10)

ans = 

whos ans

  Name      Size            Bytes  Class    Attributes

  ans       1x24                8  sym                

To convert to double numeric for use in other Matlab functions that only accept numeric types, just cast to a double

format long e
double(f21_sol)

ans = 1×24

-7.004252767019106e+01 -8.640952260733972e+01 -1.241972034940549e+02 -3.609681083944062e+02 1.391165845886463e+02 2.640676617767998e+01 -2.170360066723050e+00 -1.597473733509369e+01 -2.453099071727911e+01 -3.059052030227865e+01 -3.523107329686412e+01 -3.895083847594472e+01 -4.200308652976788e+01 -4.452284364555465e+01 -4.658211633145788e+01 -4.821601058357839e+01 -4.943481290580725e+01 -5.022750964769827e+01 -5.055781422460217e+01 -5.035013744294519e+01 -4.945652765890587e+01 -4.757812864210533e+01 -4.405409996415674e+01 -1.023378794595449e+02

[Brittany Burns]

Thank you for the explanation!

Answer 2

[Walter Roberson]

V1 = -(760*(348475767795489804396578940311480912*3^(1/2) + 168032379110159424438402050627659277))/(3*(3858879662023561551604431872957445612*3^(1/2) - 3892977404894321714789280647602584263))

V1 = -70.0425

V2 = str2sym('-(760*(348475767795489804396578940311480912*3^(1/2) + 168032379110159424438402050627659277))/(3*(3858879662023561551604431872957445612*3^(1/2) - 3892977404894321714789280647602584263))')

V2 = 

V2s = simplify(V2)

V2s = 

What are you expecting the value to simplify to? Are you looking to get the denominator in rational form?

[N,D] = numden(V2);
D1 = children(D,1);
D2 = children(D,2);
Dc = D1 - D2;
expand(N*Dc) / expand(D*Dc)

ans = 

My suspicion is that you are misunderstanding what it means to simplify an expression. Simplification of an expression does not change the value of an expression (other than it might end up removing some discontinuities.) The value of the expression involves the irrational number sqrt(3) and so is itself an irrational number -- something that is infinitely long in its decimal expansion.

I suspect what you are looking for is a decimal approximation. Which you can get using vpa() or double(), but not using simplify(). Because a decimal approximation is a different number that is only close to the original number, whereas simplify() preserves the original number.

[Walter Roberson]

dbtype foot_forces.m

1     %Finding Forces in Foot using CSA method
2     %% Set up
3     clear all; clc
4     load('scuba_motion.mat','theta*')
5     load('scuba_muscles.mat','insert*','origin*')
6     
7     % Set Constants
8     theta = theta_ankle;
9     gastro_theta = linspace(-20,50,length(theta_ankle));
10    tibialis_anterior_theta = linspace(15,20,length(theta_ankle));
11    tibialis_posterior_theta = linspace(5,15,length(theta_ankle));
12    extensor_digitorium_longus_theta = linspace(15,20,length(theta_ankle));
13    extensor_hallucis_longus_theta = linspace(15,20,length(theta_ankle));
14    peroneus_longus_theta = linspace(30,50,length(theta_ankle));
15    peroneus_brevis_theta = peroneus_longus_theta;
16    peroneus_tertius_theta = linspace(15,20,length(theta_ankle));
17    
18    %weight and height
19    fin_wt = 1.0*9.81; %fin weight (N)
20    wt = 147.5*0.453592; %total wt kg
21    ht = 64.6*2.54; %total height in cm
22    ft_wt = wt*0.0133*9.81; %foot weight (N)
23    w = fin_wt+wt; %fin plus foot weight (N)
24    
25    %cross sectional areas (CSA)
26    csa_Gastrocnemius = 22.5;
27    csa_Soleus = 25.0;
28    csa_Tibialis_Anterior = 12.5;
29    csa_Tibialis_Posterior = 12.0;
30    csa_Extensor_Digitorum_Longus = 8.0;
31    csa_Extensor_Hallucis_Longus = 6.5;
32    csa_Peroneus_Longus = 12.0;
33    csa_Peroneus_Brevis = 8.5;
34    csa_Peroneus_Tertius = 4.0;
35    
36    csa_f21 = csa_Tibialis_Anterior+csa_Extensor_Hallucis_Longus;
37    csa_f18 = csa_Gastrocnemius+csa_Soleus;
38    csa_f22 = csa_Tibialis_Posterior;
39    csa_f23 = csa_Extensor_Digitorum_Longus+csa_Peroneus_Tertius;
40    csa_f28 = csa_Peroneus_Brevis+csa_Peroneus_Longus;
41    
42    % k's
43    k21=csa_f18/csa_f21;
44    k31=csa_f22/csa_f21;
45    k41=csa_f23/csa_f21;
46    k51=csa_f28/csa_f21;
47    
48    % distances
49    
50    a=[ insert_foot(21,1) insert_foot(18,1) insert_foot(22,1) insert_foot(23,1)...
51        insert_foot(28,1) 4.5];
52    
53    %% Equations
54    
55    %Force and Moment
56    syms f21 f18 f22 f23 f28 fj
57    
58    for i=1:length(theta)
59    
60        e1 = f21.*cosd(tibialis_anterior_theta(i))+f18.*cosd(gastro_theta(i))+...
61            f22.*cosd(tibialis_posterior_theta(i))+f23.*cosd(extensor_digitorium_longus_theta(i))+...
62            f28.*cosd(peroneus_brevis_theta(i))-fj.*cosd(theta(i)) == 0; %Force in x-dir
63    
64        e2 = f21.*sind(tibialis_anterior_theta(i))+f18.*sind(gastro_theta(i))+...
65            f22.*sind(tibialis_posterior_theta(i))+f23.*sind(extensor_digitorium_longus_theta(i))+...
66            f28.*sind(peroneus_brevis_theta(i))-fj.*sind(theta(i))-w == 0; %force in y-dir
67    
68        e3 = a(1)*f21+a(2)*f18+a(3)*f22+a(4)*f23+a(5)*f28-a(6)*w == 0; %moment about ankle joint
69    
70    
71        % Forces relative to each other
72        e4 = k21*f21-f18 == 0;
73        e5 = k31*f21-f22 == 0;
74        e6 = k41*f21 -f23 == 0;
75    
76        % Solve system of nonlinear equations
77        solutions(i,:) = solve([e1, e2, e3, e4, e5, e6], [f21, f18, f22, f23, f28, fj]);
78        f21_sol(i)=solutions(i).f21;
79        f18_sol(i)=solutions(i).f18;
80        f22_sol(i)=solutions(i).f22;
81        f23_sol(i)=solutions(i).f23;
82        fj_sol(i)=solutions(i).fj;
83    
84    end

You are using lots of floating point numbers, and trig functions applied to floating point numbers, but you are also using solve().

This is a "category error". You are using the wrong function.

Floating point numbers are inherently not precise, and cannot be made precise.

Double precision point 9.81 does not represent exactly 981 / 100: it represents the (infinite) set of numbers that are between 9.81 * (1-eps/2) and 9.81*(1+eps/2) where eps is . The representative number for that set is 9.8100000000000004973799150320701301097869873046875 exactly.

Likewise 2.54 represents a set of numbers, and the representative number for that set is 2.54000000000000003552713678800500929355621337890625

If you were to multiply 9.81 * 2.54 you would not get (981 * 254) / (100 * 100) = 249174 / 10000 -- you would get 24.9174000000000006593836587853729724884033203125 exactly.

So when you make computations with floating point numbers that have even one non-zero decimal place, you are inherently saying, "I have here a number that I do not know the exact value of, but here is the value to about 16 decimal place.". And when you use 9.81 for gravity, you would be lying about it being to about 16 decimal places, as the actual force of gravity is closer to 9.80665

Your input constants are numeric garbage compared to the hypothesis that they are "exact" values.

But at the same time you are using solve(). And solve() is a function reserved for trying to find exact solutions as much as possible. So if the mathematical result is then the value from solve() is going to be 2^(1/2) not 1.4 or 1.41 or 1.4142135623730951454746218587388284504413604736328125

If your values are not precise then using solve() is like using an atomic force microscope to paint lipstick on something that is at best mostly in the porcine family, while being actively worried about mistakes caused by quantum indeterminacy. solve() is the ultimate in Garbage In Garbage Out .

If your inputs are crude numeric approximations, either just use floating point calculations, or else use vpasolve()... and later vpa() down to an appropriate number of digits. (Which, in your case, considering how low the resolution is of your coefficients, is probably at best the 10's place, if you are lucky.)