# Matlab limitation in fsolve using function input

3 views (last 30 days)
Meva on 22 Aug 2016
Commented: Meva on 31 Aug 2016
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)];

Walter Roberson on 22 Aug 2016
Edited: Walter Roberson on 22 Aug 2016
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 on 31 Aug 2016
A sincere apologise Walter. Sorry.

Alan Weiss on 22 Aug 2016
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 on 22 Aug 2016
You are right. However, it does give me again incorrect plot as attached. It does not understand my loop!