Matlab limitation in fsolve using function input

Hello,

I tried to loop for time value (T) inside my fsolve, but fsolve is pretty unforgiving.

The time loop does not seem working.

When I plot, it gives the same values (h=x(1) and theta=x(2) does not change over time which should change)!

Please see the the script that uses for loop for time (T). T is input for fsolve. :

x0 = [.1, .1];
options = optimoptions('fsolve','Display','iter');
dt=0.01;
Nt=1/dt+1;
Tarray = [0:dt:1];
T = 0;
for nt=1:Nt   
[x,fval] = fsolve(@torder1,x0,options,T)
 T=T+dt;
h(nt)=x(1);
theta(nt) = x(2);

plot(Tarray,h,'*')
hold on
plot(Tarray,theta,'+')
end

and the function for fsolve:

function F=torder1(x,T)

     x_1=[0:0.01:1];
     b=0.6;

%$ sol(1) = h; sol(2) =theta;
clear x_1;   
syms x_1 h theta kappa
f_1(x_1,h,theta,kappa) = 1/2*(1-( (h+(1-b)*theta)^2/(h+(x_1-b)*theta-kappa*x_1*(1-x_1))^2 ));
f_2(x_1,h,theta,kappa) = (x_1-b)/2*( 1-( (h+(1-b)*theta)^2/(h+(x_1-b)*theta-kappa*x_1*(1-x_1))^2 ));
kappa =1;

f_11 = 1-( (h+(x_1-b)*theta)^2/(h+(x_1-b)*theta-1*x_1*(1-x_1))^2 );
f_21 = (x_1-b)/2*( 1-( (h+(1-b)*theta)^2/(h+(x_1-b)*theta-x_1*(1-x_1))^2 ));

fint_1 = int(f_11, x_1);
fint_2 = int(f_21, x_1);

x_1=1;
upper_1=subs(fint_1);
upper_2=subs(fint_2);
clear x_1;
x_1=0;
lower_1=subs(fint_1);
lower_2=subs(fint_2);
clear x_1;
integral_result_1old=upper_1-lower_1;
integral_result_2old=upper_2-lower_2;

h0 = kappa *b*(1-b);
theta0 = kappa*(1-2*b);

integral_result_1 = subs(integral_result_1old, {h, theta}, {x(1), x(2)});
integral_result_2 = subs(integral_result_2old, {h, theta}, {x(1), x(2)});

F = [double(x(1) - integral_result_1*T^2 -h0);
    double(x(2) - integral_result_2*T^2 - theta0)];

Answer 1

[Walter Roberson]

fsolve() is for real values, but solutions to your equations are strictly complex, except at T = 0.

You might be getting false results from fsolve(), with it either giving up or finding something that appears to come out within constraints.

For example, if you

fsolve(@(x) x^2+1, rand)

then you will get a small negative real-valued answer that MATLAB finds to be within the tolerances, when the right answer should be something close to sqrt(-1)

[Meva]

Thanks Walter. Can I just omit the imaginary parts then continue looping then??

If it might give false results, how should I find a solution which method should I use? Thanks..

[Meva]

Is there any way more efficient than fsolve for complex cases like this one ? Thanks.

[Torsten]

Solve for real and imaginary part of your solution separately.

E.g. if you want to solve

z^2+1=0,

solve

real(z^2+1) = x^2-y^2+1 = 0
imag(z^2+1) = 2*x*y = 0

for

z = x+1I*y

Best wishes

Torsten.

[Meva]

Dear Torsten, I do not understand what you meant. I should repeat my query :

Is there any way for my system instead of FSOLVE?

Thanks ..

[Torsten]

No, do it in FSOLVE the way I showed you above.

Best wishes

Torsten.

[Meva]

Ok. I will separate real and imaginary parts and then I will reunify them in a loop. Thanks.

[Walter Roberson]

You do not need to separate the real and imaginary parts. You can probably use vpasolve()

X = sym('x', [1,2]);
T = 1/100;
F = torder1(X,T);
sols = vpasolve(F, X);
sols.x1
sols.x2

Note: this will find at most one solution per vpasolve, and the solutions for adjacent T values will not necessarily be "close" to each other. For example if the function is more or less symmetric around 0 in one of the variables, then you might get either root for any one fsolve().

[Meva]

Dear Walter, I do appreciate your answer. However, if I tried your answer above in the command window, I get the following error:

DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use VPA.

It does angry with double expressions in F.

But these double expressions are necessary to convert F symbolic function to a symbolic expression.

[Meva]

If I do omit double expression inside F function, and then type :

X = sym('x', [1,2]);
T = 1/100;
F = torder1(X,T);
sols = vpasolve(F, X);
sols.x1
sols.x2

This gives the following error :

CAT arguments dimensions are not consistent.

Alternatively, If I omit double expression inside F function, and then type :

X = sym('x', [1,2]);
T = 1/100;
F = torder1(X,T);
sols = vpa(F);
sols.x1
sols.x2

It still gives the CAT error as above.

So I conclude that there is a conflict between prescribing your function arguments in the type of double, and VPA and/or vpasolve.

[Walter Roberson]

I sleep. I have drive failures. I have network failures. I have appointments.

[Meva]

A sincere apologise Walter. Sorry.

Answer 2

[Alan Weiss]

I think that you need to replace your line

 [x,fval] = fsolve(@torder1,x0,options,T)

with

 [x,fval] = fsolve(@(x)torder1(x,T),x0,options)

Alan Weiss

MATLAB mathematical toolbox documentation

[Meva]

Nope!

[x,fval] = fsolve(@torder1(x,T),x0,options)

is invalid MATLAB syntex.

[John D'Errico]

But that is not what Alan suggested. There is a difference between these lines:

[x,fval] = fsolve(@torder1(x,T),x0,options)
[x,fval] = fsolve(@(x)torder1(x,T),x0,options)

Alan suggested the second, but you then tried the first.

[Meva]

You are right. However, it does give me again incorrect plot as attached. It does not understand my loop!