why is my ea not being recognized?

clc; clear all;
L = 1.25; E = 50000; I = 30000; w0 = 2.5;
syms x
f(x) = (w0/(120*E*I*L))*(-x^5+2*(L^2)*(x^3)-(L^4)*x);
df = diff(f,x); 
true = df(2); 
roots = solve(f==0, x);
true_roots = double(roots);
xl = 200;
xu = 300;
xr = (xl + xu) / 2;
es = 0.001;
xr_new(1) = xr;
if f(xl) * f(xu) < 0
    disp('The function changes sign within this bound')
    for i = 1:100
        val(i) = f(xl) * f(xr_new(i));
        if val(i) < 0
            disp('The root lies in the lower interval')
            xu = xr_new(i);
        elseif val(i) > 0
            disp('The root lies in the upper interval')
            xl = xr_new(i);
        else val = 0;
            fprintf('The approximate root of the function is %1.2f \n',xr_new(i))
            break
        end
        xr_new(i+1) = (xl + xu) / 2;
        ea = abs((xr_new(i+1) - xr_new(i))./xr_new(i+1));
        if ea < es
            break
        end
    end
else
    disp('The function does not change sign within this bound')
end

The function does not change sign within this bound

% Print the results
ea = ea*100;

Unrecognized function or variable 'ea'.

fprintf('The absolute approximate error is: %1.5f \n',ea)
n_iterations = i+1;
fprintf('The total number of iterations are: %d \n',n_iterations)
root = xr_new(1,end);
fprintf('The approximate root of the function is: %1.5f \n', root)
x1 = -1000 : 1000;
y1 = f(x1);
y = double(y1);
xmin = -1000; xmax = 1000; ymin = -0.5; ymax = 0.5;
plot(x1, y1)
hold on
plot(true_root,0,'mo','MarkerFaceColor','m')
set(gca,'XAxisLocation','origin','YAxisLocation','origin','XMinorTick','on')
xlabel('x \rightarrow')
ylabel('\uparrow f(x)')
title('Graphical root')

Answer 1

[Sulaymon Eshkabilov]

There are a few inconsistencies.

(1) The simulation in [for ... end] loop never reached due to if condition f(xl)*f(xu) <0. The actual values of f(xl)*f(xu) were in sym and greater than 0. Therefore, by changing it to double and the if condition (f(xl)*f(xu)>), you will get the simulation results.

(2) The found true roots were in sym.

(3) Similarly, y1 = f(x1) was in sym that need to be numerical. Thus, double() is needed.

(4) if val = 0 is not correct that must be val == 0.

(5) It is better (must be) to use fplot() instead of plot to plot actual faction. By this way tru solutions can be shown explicitly - see figure (2).

Here is the corrected code:

clc; clearvars;
L = 1.25; E = 50000; I = 30000; w0 = 2.5;
syms f(x)
f(x) = (w0/(120*E*I*L))*(-x^5+2*(L^2)*(x^3)-(L^4)*x);
df = diff(f,x); 
true = double(df(2)); 
roots = vpasolve(f==0, x);
true_roots =double(roots(:));
xl = 200;
xu = 300;
xr = (xl + xu) / 2;
es = 0.001;
xr_new(1) = xr;
if double(f(xl)) * double(f(xu)) >0
    disp('The function changes sign within this bound')
    for i = 1:100
        val(i) = double(f(xl)) * double(f(xr_new(i)));
        if val(i) < 0
            disp('The root lies in the lower interval')
            xu = xr_new(i);
        elseif val(i) > 0
            disp('The root lies in the upper interval')
            xl = xr_new(i);
        else val == 0;
            fprintf('The approximate root of the function is %1.2f \n',xr_new(i))
            break
        end
        xr_new(i+1) = (xl + xu) / 2;
        ea = abs((xr_new(i+1) - xr_new(i))./xr_new(i+1));
        if ea < es
            break
        end
    end
else
    disp('The function does not change sign within this bound')
end

The function changes sign within this bound

The root lies in the upper interval
The root lies in the upper interval
The root lies in the upper interval
The root lies in the upper interval
The root lies in the upper interval
The root lies in the upper interval
The root lies in the upper interval
The root lies in the upper interval

% Print the results
ea = ea*100;
fprintf('The absolute approximate error is: %1.5f \n',ea)

The absolute approximate error is: 0.06515 

n_iterations = i+1;
fprintf('The total number of iterations are: %d \n',n_iterations)

The total number of iterations are: 9 

root = xr_new(1,end);
fprintf('The approximate root of the function is: %1.5f \n', root)

The approximate root of the function is: 299.80469 

x1 = -1000 : 1000;
y1 = double(f(x1));
xmin = -1000; xmax = 1000; ymin = -0.5; ymax = 0.5;
figure(1)
plot(x1, y1)
hold on
plot(true_roots,0,'mo','MarkerFaceColor','m')
set(gca,'XAxisLocation','origin','YAxisLocation','origin','XMinorTick','on')
xlabel('x \rightarrow')
ylabel('\uparrow f(x)')
title('Graphical root')
hold off

figure(2)
Fun = @(x) (w0/(120*E*I*L))*(-x.^5+2*(L^2)*(x.^3)-(L^4)*x);
fplot(Fun, [-1.5, 1.5], 'kx-')
hold on
plot(true_roots,0,'mo','MarkerFaceColor','m')
set(gca,'XAxisLocation','origin','YAxisLocation','origin','XMinorTick','on')
xlabel('x \rightarrow')
ylabel('\uparrow f(x)')
title('Graphical root')

[Clinton Mitchell]

Thank you so much

[Sulaymon Eshkabilov]

Most welcome!

[Walter Roberson]

The simulation in [for ... end] loop never reached due to if condition f(xl)*f(xu) <0. The actual values of f(xl)*f(xu) were in sym and greater than 0. Therefore, by changing it to double and the if condition, you will get the simulation results.

No, this is incorrect logic. f(xl) and f(xu) were both negative, so their product was positive. A positive product of the values indicates that the two are the same sign, which means there are an even number of sign changes between the endpoints (possibly zero sign changes), which is something that the bisection algorithm cannot deal with.

You "fixed" the problem by changing the test for f(xl)*f(xu) < 0 [a correct test] to f(xl)*f(xu) > 0 [an incorrect test]. The real problem was that the initial xl and xu were way outside of the places where the sign changed.

(2) The found true roots were in sym. Therefore vpasolve() is needed.

This is incorrect as well.

solve() looks for analytic roots. For polynomials it always generates exact solutions or placeholder structures standing in for exact solutions. For non-polynomials that involve exp() or trig functions, it tries to determine whether it can find all solutions analytically, and if so does so. For non-polynomials that are not one of the forms it can deal with analyatically, and where the number of equations matches the number of variables to be solved for, solve() calls vpasolve() to return a single solution; if the number of equations does not match then solve returns unresolved or empty.

vpasolve() does not look for analytic roots. However, for polynomials it always generates the full list of solutions as numeric approximations. For non-polynomials, it uses a modified a Newton method to try to find a single numeric solution. If the number of equations does not match the number of variables, vpasolve() generates an error.

Both solve() and vpasolve() return solutions of class sym . vpasolve() returns symbolic floating point numbers -- which are still class sym . solve() returns analytic expressions when possible . In the cases where those expressions do not have any free variables, vpa() can be used to approximate them as symbolic floating point numbers, or double() can be used to approximate them as double precision. Both of them return sym so double() of either of them is the same effect.

solve() is required if you are solving for fewer variables than exist in the expression. vpasolve() is usually faster than solve(), potentially much faster. But there are cases where solve() is able to find solutions that vpasolve() is not able to find, in cases where continuity is an issue, or cases where the function is very steep. In cases that do not involve inequalities, if vpasolve() is able to find a solution but solve() says that solutions do not exist (not just that it cannot find an explicit solution) then the large majority of the time, the numeric solution is a false root.

Answer 2

[Voss]

ea is not defined because f(xl) * f(xu) >= 0

L = 1.25; E = 50000; I = 30000; w0 = 2.5;
syms x
f(x) = (w0/(120*E*I*L))*(-x^5+2*(L^2)*(x^3)-(L^4)*x);
xl = 200;
xu = 300;

f(xl) * f(xu)

ans = 

so the code inside "if f(xl) * f(xu) < 0" does not execute; only the code inside the "else" block executes.

Note the displayed message: "The function does not change sign within this bound"

Answer 3

[Walter Roberson]

if f(xl) * f(xu) < 0

In the case where that does not hold true, the for loop is never entered and ea is never assigned to.

Notice you get the output message "The function does not change sign within this bound"

[Walter Roberson]

the true_roots shows you that the actual zeros are contained within [-2 2], no-where near the 200 to 300 you are searching.

[Clinton Mitchell]

So how do I fix it then

[Walter Roberson]

Change

else
    disp('The function does not change sign within this bound')
end

to

else
    disp('The function does not change sign within this bound')
    ea = inf;
end