How can I adjust the tolerance in vpasolve?

I am trying to solve a (complex) system of two equations with two unknowns using vpasolve. However, there does not seem to be a exact value which satisfies both equations. Hence, I would like to obtain the closet possible value to the solution adjusting the tolerance. However, it seems like there is not this possibility.

Could you please give me hand?

Thanks

Here are my equations, which depend on t and n:

 2.991228070175400e-09=0.0049*(((n-1.32)*n*(1.209e-04)*(1-exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5)*t)))/(1+n^2*(1.209e-04)*(1-exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5).*t))+1.32^2.*(1.209e-04)*exp(-2*2*pi/(1550e-9)/2.7*(2.7-1.32)^(-0.5).*t)));
 2.357536764705900e-09=0.0049*(((n-1.32).*n.*(1.3562e-04)*(1-exp(-2*2*pi/(1550e-9)/3.25*(3.25-n)^(-0.5).*t)))./(1+n^2.*(1.3562e-04)*(1-exp(-2*2*pi/(1550e-9)/3.25_2*(3.25_2-n)^(-0.5).*t))+1.32^2.*(1.3562e-04)*exp(-2*2*pi/(1550e-9)/3.25*(3.25-1.32)^(-0.5).*t)));

[John D'Errico]

You claim that they depend only on t and n. But I see B_2 in there.

[pepe juan]

Sorry, forgot to remove it. I have just corrected it.

[John D'Errico]

Looks like Stephen removed the probably spurious _2 characters in there.

Answer 1

[John D'Errico]

Anyway, even ignoring the fact that there is a variable B_2 in there that has never been defined...

I see too often people thinking that a tighter tolerance will help them. Just demanding a finer solution does not always help. For example, consider just one of the numbers in those expressions. You may think that the symbolic toolbox understands what you have typed. But did it really?

Look at the number 2.991228070175400e-09. sym turns it into a rational fraction.

N = sym(2.991228070175400e-09)
N =
3616172846391081/1208925819614629174706176
>> vpa(N,60)
ans =
0.00000000299122807017540003409435527014816058599677717211307026445866

But as you see, when we look at that fraction, re-expanded to a decimal, we see that it was not the original number you provided. In fact, the fraction approximates the number that is stored in the floating point number. Here is that floating point number, written out to many digits.

sprintf('%1.100f',2.991228070175400e-09)
ans =
0.0000000029912280701754000340943552701481605859967771721130702644586563110351562500000000000000000000

In both cases, MATLAB has actually represented the number that you typed with 16 significant digits or so, but stored it in binary form, so not exactly. The number you typed in is not the number that was stored. The point is, then you write a number like that, MATLAB first represents it as a double precision number, so stored in binary form, using only 53 binary bits of precision. So about the first 16 significant decimal digits will be correct. But in fact, you cannot represent a number like that decimal exactly in binary.

So while you THINK you provided some long set of equations that you want an exact solution to. In fact, you provided a set of equations with only very approximate coefficients. Those coefficients are not represented exactly, because they were provided in double precision form.

In the end, you are trying to find an extremely precise solution to the wrong system of equations. Feel free to try though. You are probably just wasting a lot of CPU cycles. This may help you:

help digits

But is that REALLY the number you wanted? Where exactly did you get the number 2.991228070175400e-09? In fact, it was surely computed from something. Was that computation exact? I think you are kidding yourself. Now consider the same questions for EVERY floating point number that you wrote in those long equations. Arbitrarily, take another number in there:

N = sym(6.045e-05)
N =
8920845434045939/147573952589676412928

vpa(N+N,100)
ans =
0.000120899999999999997811299390360062488980474881827831268310546875

As you see, what you wrote is not how those equations will be evaluated.

[pepe juan]

Firstly, many thanks for your in-depth and constructive comment.

Yes, the number 2.991228070175400e-09 was computed from a fitted function (here another approximation).

But as you point out, I do not need that much precision. In fact a solution which satisfies the equations so that I obtain 2.99e-09 would be enough for what I am aiming to do.

[John D'Errico]

But cranking down in the tolerance will not help it to solve an apparently intractable problem.

I'd generally suggest first plotting things. WHENEVER you have a problem like this, ok, even before you have a problem like this, PLOT IT!!!!!!!!! Never just throw some problem into a solver without understanding what you are trying to solve.

For example, if I take the first equation, and convert it into something that MATLAB can use in a simple way...

E1 = @(n,t) 0.0049.*(((n-1.32).*n.*(6.045e-05+6.045e-05)*(1-exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5).*t)))./(1+n^2.*(6.045e-05+6.045e-05)*(1-exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5).*t))+1.32^2.*(6.045e-05+6.045e-05)*exp(-2*2*pi/(1550e-9)/2.7*(2.7-1.32)^(-0.5).*t))) - 2.991228070175400e-09;

Then I would see that the problem appears to go through a singularity around n = 2.5 or so. Inspecting the equation, on would see that above n=2.7, there is a problem, where the solutions will suddenly turn into complex numbers. So n has a maximum of 2.7. As well, if you recreate the above plot, carefully rotating it, then you will see there appear to be solutions, where the E1 surface does cross the axis at zero. But it does so just barely, and the surface appears to hardly be a function of t at all.

In fact, there was no need to do anything symbolic.

But, when I try to do the same thing with the second equation, I see a multitude of errors. There are lots of places where MATLAB fails, trying to evaluate sub-expressions like: 3.25_2

So my guess is this is not correct code. Just a guess. IF I were to guess, then I MIGHT decide that maybe you just are not good at editing your code, and that maybe you replaced a variable name with some constant, like 3.25. but that the variable name ended with "_2". So, if I delete the two instances where that happens...

E2 = @(n,t) 0.0049*(((n-1.32).*n.*(1.3562e-04)*(1-exp(-2*2*pi/(1550e-9)/3.25*(3.25-n)^(-0.5).*t)))./(1+n^2.*(1.3562e-04)*(1-exp(-2*2*pi/(1550e-9)/3.25*(3.25-n)^(-0.5).*t))+1.32^2.*(1.3562e-04)*exp(-2*2*pi/(1550e-9)/3.25*(3.25-1.32)^(-0.5).*t))) - 2.357536764705900e-09;

So, lets change the axes, replot, then change the view.

ezsurf(E1,[-5,2.7,-10,10])
hold on
ezsurf(E2,[-5,2.7,-10,10])
view(0,0)

So, they are clearly two subtly different surfaces. Neither of which appears to be a function of t to any significant amount. The two surfaces both cross Z==0 along two nearly linear paths, both close to zero.

So all of this was done without the benefit of any symbolic tools.

So, lets look more closely yet, expanding on the region of interest in t, but narrower yet in n.

ezsurf(E1,[-.15,1.4,-1000,1000])
hold on
ezsurf(E2,[-.15,1.4,-1000,1000])
view(0,0)

Again, it appears, even allowing t to vary widely, that the two surfaces are virtually independent of t. By viewing the surface head on, we see a straight line only. So t impacts virtually nothing in the shape of that surface.

[John D'Errico]

The previous comment started getting too long to edit, so I'll continue with another comment. The point is, you need to look next at how n and t enter the expression.

So in E1, I see t only inside the exponentials, where it becomes obvious that only when t is rather large, on the order of 1e10, will anything happen at all as a function of t.

exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5).*t)

Think about how that expression will be evaluated, and what it will take for that expression to be significant as a function of t.

Going back to symbolic form, we can see

E1_Sym = 0.0049.*(((n-1.32).*n.*(6.045e-05+6.045e-05)*(1-exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5).*t)))./(1+n^2.*(6.045e-05+6.045e-05)*(1-exp(-2*2*pi/(1550e-9)/2.7*(2.7-n)^(-0.5).*t))+1.32^2.*(6.045e-05+6.045e-05)*exp(-2*2*pi/(1550e-9)/2.7*(2.7-1.32)^(-0.5).*t)))-2.991228070175400e-09;

E2_Sym = 0.0049*(((n-1.32).*n.*(1.3562e-04)*(1-exp(-2*2*pi/(1550e-9)/3.25*(3.25-n)^(-0.5).*t)))./(1+n^2.*(1.3562e-04)*(1-exp(-2*2*pi/(1550e-9)/3.25*(3.25-n)^(-0.5).*t))+1.32^2.*(1.3562e-04)*exp(-2*2*pi/(1550e-9)/3.25*(3.25-1.32)^(-0.5).*t))) - 2.357536764705900e-09;

vpa(subs(E1_Sym,n,1.32))
ans =
-0.00000000299122807017540003409435527014816058599677717211307026445865631103515625

So, when n=1.32 or 0, this mess is actually an exact constant, independent of t. For other values of n however, there is that tiny exponential aspect of t.

But, lets try fixing the value of t, at say, 1. Then solve for n.

vpasolve(subs(E1_Sym,t,1),n)
ans =
-0.00381417066708864375383004372895748541442098531057225352112150250677438336862087403638915176878066005

vpasolve(subs(E2_Sym,t,1),n)
ans =
-0.002682150613590733634837035120794878802053116920717520146875039220026227663376508594850127789943593552

So if I fix t at 1, then I get different solutions for n from the two equations.

So, where does t become significant? It turns out those exponential forms become significant only near t==0, for positive values of t.

So lets go in really fine, near one of the near solutions around n=0.

ezsurf(E1,[-.1,0.1,0,1e-6])
hold on
ezsurf(E2,[-.1,0.1,0,1e-6])

What we are looking for here are the places where the lines in those overlapping surfaces are fuzzy, and darker than otherwise, AND where the surfaces cross z==0. What is where your solution must lie. IF any solution exists.

So I'll now try using vpasolve, but with some intelligenece applied to the solution, using good starting values for n and t.

nt = vpasolve(E1_Sym,E2_Sym,[n,t],[-1.e-5,1.e-7])
nt = 
    n: [0x1 sym]
    t: [0x1 sym]

This suggests the surfaces may possibly never touch down at zero in the same place near n==0.

nt = vpasolve(E1_Sym,E2_Sym,[n,t],[1.33,1.e-7])
nt = 
    n: [1x1 sym]
    t: [1x1 sym]

nt.n
ans =
1.39114688409837609064028449811

nt.t
ans =
0.0000000199546731461591267421429033937

So it did find a solution near n=1.32. Was it successful?

subs(E1_Sym,[n,t],[nt.n,nt.t])
ans =
8.75811540203010669327330989556e-47

subs(E2_Sym,[n,t],[nt.n,nt.t])
ans =
2.62743462060903200798199296867e-46

So via a rather long, arduous process, there appears to be a solution. Nothing more than common sense was needed, in choosing intelligent starting values for the solver, and in understanding the fundamental shapes of the surfaces involved. No cranking down on the tolerances though. I did decide to use vpasolve because the numbers were so nasty, and because it would be so easy to have problems in double precision arithmetic. I might possibly have used fsolve with some careful scaling of the problem.

[John D'Errico]

In the end, you need to next do some sensitivity analysis of that solution. You might want to decide if subtle perturbations in the coefficients of that function are able to then induce significant changes in the solution. How well did you know those coefficients?

I would NEVER trust that the above solution is really a solution until I knew enough about the problem, what a solution means, where the equations came from, and how well I know that mess of coefficients.

Answer 2

[Touseef]

Each question output with command fprintf and decimal places of solution should be one greater than

the tolerance level. (For example solution will be displayed 0.7391 for accuracy of 10􀀀3)???

[John D'Errico]

Please, don't ask a question as an answer to another question. Anyway, this question does not even make sense. fprintf cannot possibly know what was the tolerance employed to result in a given number. As far as fprintf knows, it is just a number at the point where fprintf is used.