# float point arithmetic and vpasolve

6 views (last 30 days)
Meva on 30 Aug 2016
Commented: John D'Errico on 30 Aug 2016
Hi,
I have symbolic integrations. (int_1 , int_2 inside t.m).
I need to solve an equation system with these integrations. First, I used fsolve. But it did not work since the integrations return complex values. Then I used vpasolve. The question is How trustable is vpasolve? If I use tlast =0.7, Matlab says :
Error using plot Vectors must be the same length.
If I use tlast= 0.77 it does not give me an error. However, the plot has strange looking.
I suspect that it does not calculate F when T=0.7 as the plot for harray and thetaarray versus Tarray should be smooth.
Is there any wayout for my problem. Codes are below
clear all; close all; clc;
x0 = [.1, .1];
options = optimoptions('fsolve','Display','iter');
dt=0.07;
M = .1;
g = .01;
Gamma = 0.2;
I = Gamma*M;
tlast = 0.7;
Nt=tlast/dt+1;
Tarray = [0:dt:tlast];
kappa = 1; b=0.6;
h0 = kappa *b*(1-b);
for nt=1:Nt
T = (nt-1)*dt;
X = sym('x', [1,2]);
F = t(X,T,M,g,Gamma);
sols = vpasolve(F, X);
h_actual=sols.x1
theta_actual=sols.x2
h = double(sols.x1)
theta = double(sols.x2)
harray(nt) = h
thetaarray(nt) = theta
end
figure(1)
subplot(2,1,1)
plot(Tarray,harray,'k')
subplot(2,1,2)
plot(Tarray,thetaarray,'k')
The function I used is
function F=t(x,T,M,g,Gamma)
x_1=[0:0.01:1]
b=0.6;
syms x_1 h theta
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 ));
int_1 = int(f_11, x_1);
int_2 = int(f_21, x_1);
x_1=1;
upper_1=subs(int_1);
upper_2=subs(int_2);
clear x_1;
x_1=0;
lower_1=subs(int_1);
lower_2=subs(int_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 = [vpa(x(1) - (1/2*integral_result_1 - M*g)*T^2 -h0); vpa(x(2) - (1/(2*Gamma)*integral_result_2)*T^2 - theta0)];
Before using vpasolve I did use fsolve. However, it did not solve my problem for complex values. The function for fsolve does not allow to take T for inputs. Even if it allows, when I plot harray versus Tarray I get wrong plot. Is MAtlab capable of solving the system I have. Thanks..
John D'Errico on 30 Aug 2016
Luckily, hammers have no feelings to be hurt. If they did, then I'd keep a careful watch on your thumbs, as a vindictive hammer would be unpleasant to have around.

John D'Errico on 30 Aug 2016
Edited: John D'Errico on 30 Aug 2016
I've tried to delete some of the extraneous fluff you have in that code, but even there, you have some serious potential problems.
dt=0.07;
M = .1;
g = .01;
Gamma = 0.2;
I = Gamma*M;
tlast = 0.7;
Nt=tlast/dt+1;
Tarray = [0:dt:tlast];
kappa = 1; b=0.6;
h0 = kappa *b*(1-b);
for nt=1:Nt
T = (nt-1)*dt;
X = sym('x', [1,2]);
F = t(X,T,M,g,Gamma);
For example, there is no reason to assign a value to x_1, if two lines later, you will re-define it as symbolic. There is no need to clear variables before you re-assign them to have new values.
Next, things like this will get you into deep trouble, because soon your next anguished question will be why floating point arithmetic does not work properly.
Nt=tlast/dt+1;
Here tlast was 0.7, and dt was 0.07. While it may indeed happen that 0.7/0.07 is exactly 10, that is not assured in floating point arithmetic. In fact, that result need not even be an integer. So, just for kicks...
r = 0.7/0.07
r =
10
It LOOKS like 10. But is it?
r == 10
ans =
0
Nope.
sprintf('%1.100f',r)
ans =
9.9999999999999982236431605997495353221893310546875000000000000000000000000000000000000000000000000000
The problem is, numbers like 0.7 are not exactly representable in floating point arithmetic. And that number is not exactly 10 times 0.07.
Welcome to the wonderful, wacky world of floating point.
So it is better for you to learn to use tools like linspace. It is always dangerous to do things like
0:.07:0.7
So,for example, what happens if we do this?
0:1:r
ans =
0 1 2 3 4 5 6 7 8 9
BUT r is 10. I think!
r
r =
10
Again, r was not truly 10.
Far safer and better is to write it as:
Nt = 11;
Tarray = linspace(0,tlast,Nt);
dt = tlast/Nt;
linspace will always get it right. (In fact, this may be part of the reason why you had troubles.)
Next, I see places where you write this:
int_1 = int(f_11, x_1);
int_2 = int(f_21, x_1);
x_1=1;
upper_1=subs(int_1);
upper_2=subs(int_2);
Using subs like that can be dangerous, because subs might grab all sorts of strange stuff to substitute. After all, I've already seen that you have a tendency to define variables multiple times. You might not realize what you put in there.
Far safer is to use subs like this:
upper_1=subs(int_1,x_1,1);
upper_2=subs(int_2,x_1,1);
This will explicitly force the desired substitution, instead of having subs go searching through your workspace for whatever it might find.
The point is, you have lots of stuff in there that I would need to check if I were you. Does it do exactly what you think it did? What you wanted it to do? When you write code, check every single line. At least do that until you are confident in your coding skills, and even then it is still a good idea.
vpasolve does work correctly, IF you use it correctly.

### Categories

Find more on Solver Outputs and Iterative Display in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!