How to avoid complex numbers out of Newton-Raphson?

I have to calculate the displacement of a beam and I get complex numbers as a solution which is a bit unrealistic.

The code is supposed to converge to the values of horizontale_displacement and vertical_displacement. Therefore over the Newton Raphson iteration the values of k and psi are getting changed for reaching that. Does anybody have a solution how I can achieve this? I get the complex numbers out of the calculations for fxpo. Also k can only reach values 0<k<1, but with right now I am at -4 +0.04i.

I am quite new to Matlab so sorry if the code is a bit messy.

clear
clc
format short g
syms R Phi1 s k psi
% Beam Characteristics
L = 5e-3;
w = 5e-6;
t = 25e-6;
I = w*t^3/12;
E = 169E9;
% Aimed displacement for calculations
vertical_displacement =200e-6;
horizontale_displacement =10e-6;
% calculation of the angles phi1 and phi2
p1 = solve(sin(Phi1) == 1/k*cos(psi/2), Phi1);
phi1 = p1(1);
phi2 = pi - phi1;
%Elliptical Integration
m = k^2;
f1 = ellipticF(phi1,m);
f2 = ellipticF(phi2,m);
F1 = vpa(f1,15);
F2 = vpa(f2,15);
e1 = ellipticE(phi1,m);
e2 = ellipticE(phi2,m);
E1 = vpa(e1,15);
E2 = vpa(e2,15);
alpha = (F2 - F1)^2;
% Displacement due to bending
bb =@(k, psi) -1/sqrt(alpha)*(sin(psi)*(2*E2 - 2*E1 - 2*F2 + 2*F1) + 2*k*cos(psi)*(cos(phi1) - cos(phi2)))*L;
ab =@(k, psi) -1/sqrt(alpha)*(cos(psi)*(2*E2 - 2*E1 - 2*F2 + 2*F1) + 2*k*cos(psi)*(cos(phi1) - cos(phi2)))*L;
%Displacement due to axis stress
Th = 2* (asin(k*sin(phi2) - asin(k*sin(phi1))));
s = 1/sqrt(alpha)*(F2-F1)*L;
lambda = sqrt(w*t*L^2/I);
fa=  cos(psi-Th)*cos(Th)*(s/L);
fb=  cos(psi-Th)*sin(Th)*(s/L);
int_fa = int(fa,0,1);
int_fb = int(fb,0,1);
ba =@(k,psi) alpha/lambda^2 *L * int_fb;
aa =@(k,psi) alpha/lambda^2 *L * int_fa;

% Formula for Displacements
Vertical_def = ab(k,psi) + aa(k,psi);
Horizontal_def = bb(k,psi) + ba(k,psi);
f = Vertical_def;
g = Horizontal_def;
% Jacobi
df = [diff(f, k) diff(f,psi)];
dg = [diff(g, k) diff(g, psi)];

% initial guesses for k and psi
k0=0.1;
psi0 =pi/5;
w = [k0 psi0];
v = w';
tol = 1e-6;
maxiter = 50;
k = v(1);
psi = v(2);
disp('')

% Displacements with initial guesses for k and psi
fun_v1= vpa(subs(Vertical_def, {'k' 'psi'},{k psi}),15);
fun_v2= vpa(subs(Horizontal_def, {'k' 'psi'},{k psi}),15);
fun_xo = [fun_v1; fun_v2];
disp('Total vertical displacement: ')
disp(fun_v1)
disp('Total horizontal displacement: ')
disp(fun_v2)

fun_1pv = vpa(subs(df, {'k' 'psi'},{k psi}),15);
fun_2pv = vpa(subs(dg, {'k' 'psi'},{k psi}),15);
fxpo = [fun_1pv; fun_2pv];
disp('Jacobi (1,1) ')
disp(fun_1pv(1))
disp('Jacobi (1,2) ')
disp(fun_1pv(2))
disp('Jacobi (2,1) ')
disp(fun_2pv(1))
disp('Jacobi (2,2) ')
disp(fun_2pv(2))
%Newton-Raphson
if det((norm(fxpo))) ==0
    disp('Error: determinatn of Jacobian is zero, try another initial point') 
else
    disp('check')
    iter = 0;
    while ((vpa(norm(fun_xo), 15))>tol) && (det(norm(fxpo))~=0) && (iter<maxiter)&& (fun_v1~=vertical_displacement)&& (fun_v1~=horizontale_displacement)
        iter = iter+1;
        vn = v - fxpo\fun_xo;
        v = vn;
        k = v(1);
        psi = v(2);   
        fv1 = vpa(subs(Vertical_def, {'k' 'psi'},{k psi}),15);
        fv2 = vpa(subs(Horizontal_def, {'k' 'psi'},{k psi}),15);
        fun_xo = [fun_v1; fun_v2];
        
        fun_1pv = vpa(subs(df, {'k' 'psi'},{k psi}),15);
        fun_2pv = vpa(subs(dg, {'k' 'psi'},{k psi}),15);
        fxpo = [fun_1pv; fun_2pv];
        if det(fxpo)==0
            disp('Error: determinatn of Jacobian is zero, try another initial point')
        end
    end
    root = v;
end

k = root(1)
psi = root(2)
fin_vert = vpa(subs(Vertical_def, {'k' 'psi'},{k psi}),15)
fin_hor = vpa(subs(Horizontal_def, {'k' 'psi'},{k psi}),15)

Thanks in advance

[Torsten]

Write out f and g and see whether these functions can produce complex results for some inputs of k and psi (e.g. x^a for a < 0, log(x) for x < 0 etc.)

I cannot check this because your code takes too long under MATLAB online.

[Dominik]

So first of all thanks a lot for the the quick response. The problem is, that coming from the diff comand over the elliptical integrals. I get quite a few of these things in the full jacobi and their result is imaginary. is there a way to keep them real?

ellipticE(asin(cos(psi/2)/k), k^2)

ellipticF(asin(cos(0.5*psi)/k), k^2)

[John D'Errico]

Simple. Don't pass in numbers that would yield complex results. For example, if you are taking a square root, and you want to keep the result real, then what would you expect IF you compute sqrt(-1)? Is there some magical way to "Keep it real"? Of course not. The only solution is to not try to evaluate the respective functions outside of the domains where they will return real results.

For example, consider the domain of asin(x). As long as x remains in the interval [-1,1], asin(x) is ALWAYS real. But if you try to evaluate asin(x) at x==2, then you will expect to see a complex result. There is simply no real value that returns a result where sin(u)=2.

asin(2)

ans = 1.5708 - 1.3170i

So the answer is simple. Avoid arguments that will yield complex results.

As far as those function fragments go, look carefully at what happens for various values of k and psi.

[Dominik]

Thanks everybody. It was a combination of some errors. I switched the vpa to evals and also adjusted my intial values as recommended. It was simpler than I thought

[Walter Roberson]

eval() is not documented for symbolic expressions or symbolic functions. There are circumstanes under which eval() of a symbolic expression will produce an error or unexpected value.

syms x
f = piecewise(-1 <= x & x <= 1, x^3, 1/sqrt(x))

f = 

x = 1/2

x = 0.5000

char(f)

ans = 'piecewise(x in Dom::Interval([-1], [1]), x^3, symtrue, 1/x^(1/2))'

eval(f)

Error using sym/eval
Invalid expression. Check for missing multiplication operator, missing or unbalanced delimiters, or other syntax error. To construct matrices, use brackets instead of parentheses.

Notice that 'x in Dom::Interval' is not valid MATLAB code, but eval() tries to execute it as-if it were valid MATLAB code.

Don't Do That.

Use subs() and possibly double() or vpa() afterwards.

Answer 1

[Walter Roberson]

Some of the issues are hidden because you are using vpa() in a number of places.

When you take a numeric approximation to a symbolic expression, because the computation has to proceed in finite precision and in some order of operations, the numeric approximation can sometimes lead to complex components that do not exist in theory. Consider for example the formal formula for roots of a cubic: two of the three roots are defined in terms of sqrt(-1) times something, and in the case where all three roots are real-valued, the complex components in theory exactly balance out to 0... but in practice they do not always exactly balance out.

I you remove the vpa() and set up your code to use rational coefficients instead of floating point numbers, so that you get indefinitely precise computations, then if you go down into the calculations, you have sub-expressions such as asin((5*2^(1/2)*(5^(1/2) + 5)^(1/2))/2) which is complex-valued, the ellipticF around that is then forced to return a complex value. That value then gets squared and just happens to balance out to exactly negative real (with no complex component)... but the negative value is then raised to ^(-1/2) and square root of negative goes complex valued again. So your code is going in and out of complex valued. I did not chase the code far enough to determine whether all of the expressions are "theoretically" real-valued after all of the manipulation being done, but by the time the initial fun_1pv is calculated, you are already dealing with expressions with complex components... it is not surprising that the situation would just get worse afterwards.