Roots of a non-linear equation using Matlab

Question

I am using Matlab to find the roots of a non-linear function. The equation is lengthy and I have used another .m to save the function, the code for which goes like

function x_c = f_x_c(s,H,VA,Lo,qc,EA,NF,Sj,Fj)

if (s < 0) || (s > Lo);
    disp('The value of s is invalid')
    disp(['s = ' num2str(s)]);
    return
end

C1 = H/qc;
if NF == 0
    n = 0;
    sn = 0;
    sum_Fj = 0;
end
if NF >= 1
    Sj_Q = [0; Sj; Lo];
    %Determine n and sn if 0 <= s < Lo:
    if s < Lo
    STOP = 0;
    k = 0;
        while STOP == 0
        k = k + 1;
            if (s >= Sj_Q(k,1)) && (s < Sj_Q((k + 1),1))
            STOP = 1;
            end
        end
    n = k - 1;
    sn = Sj_Q(k,1);
    end

    %Determine n and sn if s = Lo:
    if s == Lo
        n = NF;
        sn = Sj(NF,1);
    end

sum_Fj = sum(Fj(1:n,1));
end


x_c = (H/EA)*s;
x_c = x_c + C1*asinh((qc*s - VA + sum_Fj)/H) + ...
- C1*asinh((qc*sn - VA + sum_Fj)/H);

for j = 1:n
    sk = Sj_Q((j + 1),1);
    sk_1 = Sj_Q(j,1);
    sum_Fj = sum(Fj(1:(j - 1)));
    x_c = x_c + ...
    + C1*asinh((qc*sk - VA + sum_Fj)/H) + ...
    - C1*asinh((qc*sk_1 - VA + sum_Fj)/H);
end

The variable is H here. There is no problem with the code because it returns me that lengthy equation when I type the following in the main file.

syms x
equation = -(XB - XA) + f_x_c(s,x,VA,Lo,qc,EA,NF,Sj,Fj);  %Replaced H with variable H and all other arguments are pre-defined

Now, I want to solve this equation near H0. When I put fzero(@(x)equation, H0), it gives me an error which goes like

Undefined function 'isfinite' for input arguments of type 'sym'.

Error in fzero (line 308)
    elseif ~isfinite(fx) || ~isreal(fx)

Error in main (line 50)
fzero(@(x)equation, H0)

How can I solve this problem?

EDIT:

The equation has at least one root because if I use ezplot to plot the function, I get the following figure.

Answer 1

I see, I got to almost the same place as you, however, I got to line 309, since apparently using isfinite on a sym does not cause a crash in matlab 2014b. However, isfinite returns false, still giving an error. But to the point: it seems that fzero is not supposes to be used for symbols. Are you sure that symbols are necessary here? Looking at the function call it does not seem to be neccessary, since the output seems as if it is supposed to be a numeric anyway. Also you may argue that

syms x; fzero(@(x) x^2-1,1)

works, but then the x in @(x) have higher priority which is not a symbol (in programming we say that the variable have a different scope).

If it is important that x is a sym, you should use the equation solver solve instead that works for symbolic output.

syms x;
equation = x^2-3;
solve(equation,'x')

However, this may fail or give you incredibly long answers for complicated functions (and also for expressions not having nice fractional answers try syms x; equation = x^3-3.17+1.37*x;). It is also slower and is thus not recommended if you do not have arbitary constants in the expression (eg. f = x^2-a <=> x = +- a^(1/2), where a is to be defined later or you want to do something to the solution for multiple values of a).

Answer 2

I do not know what caused the problem as I am not familiar with the symbolic package in Matlab.

However you should be able to solve it fairly easy using some numerical solver. From what I can read of your code you want to solve the above equation equal to zero, this is exactly what Newton Rhapson is doing. You can look up the method here: http://en.wikipedia.org/wiki/Newton's_method

As you probably do not know the exact derivative of your function you can estimate it using first order approximation, where you simply use the definition of differentials but as we cannot let h go to 0, we choose h very small, in matlab I am typically using sqrt(eps). So that the approximation becomes: f'(x) = (f(x+sqrt(eps))-f(x))/sqrt(eps).

Otherwise you can use my method, which works in 1 dimension, you can find it here: http://imada.sdu.dk/~nmatt11/MM533-2014/ You have to download fpisimple.m and mynewton.m as newtons method is just an application of fix point iteration it is build on top of that.