float point arithmetic and vpasolve

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]

How trustable is vpasolve? Maybe you should first ask how trustable is the code that you wrote? I'd trust vpasolve far more than I'd trust the code of someone who has no idea how to use it.

For example, you define variables, then two lines later you define them again, differently. For some reason, you seem to think you need to clear variables, before you re-assign them to have new values.

So if I were you, I'd start reading the tutorials. Learn about variables, about the various classes for your variables, about functions.

[Stephen23]

"The question is How trustable is vpasolve?"

The question is: How trustable is a hammer ?

Like any tool, if you give a hammer to someone who knows what they are doing then it can be very useful indeed. Give the same hammer to a beginner who is not sure what they are doing or even how to hold it properly, and you will get a few broken mirrors and some nice holes in the wall...

And then they will blame the hammer! Like this:

"vpasolve does not work correctly"

There is nothing wrong with not knowing how to use a hammer and taking the time to learn how it should be used and making a few mistakes along the way (we all do this!), but blaming the hammer for the holes in the wall is not completely fair on the poor hammer...

[John D'Errico]

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.

Answer 1

[John D'Errico]

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.