MATLAB Answers

Matlab limitation in fsolve using function input

14 views (last 30 days)
Meva
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)];

  0 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
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)

  11 Comments

Show 8 older comments
Meva
Meva on 25 Aug 2016
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
Walter Roberson on 27 Aug 2016
I sleep. I have drive failures. I have network failures. I have appointments.

Sign in to comment.

More Answers (1)

Alan Weiss
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

  3 Comments

Meva
Meva on 22 Aug 2016
Nope!
[x,fval] = fsolve(@torder1(x,T),x0,options)
is invalid MATLAB syntex.
John D'Errico
John D'Errico on 22 Aug 2016
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
Meva on 22 Aug 2016
You are right. However, it does give me again incorrect plot as attached. It does not understand my loop!

Sign in to comment.

Sign in to answer this question.


Translated by