I want to solve this equation. please help me

3.8 = a (1 - exp(b×8.2)) + c (exp(d×8.1) - 1);

27.9 = a (1 exp(b×63)) + c (exp(d×35.667) - 1);

23.5 = a (1 - exp(b×59.5)) + c (exp(d×59.5) - 1);

13.7 = a (1 - exp(b×35.7)) + c (exp(d×35.7) - 1)};

Answer 1

[madhan ravi]

syms a b c d

eq1 = 3.8 == a*(1 - exp(b*8.2)) + c*(exp(d*8.1) - 1);

eq2 = 27.9 == a*(1 - exp(b*63)) + c*(exp(d*35.667) - 1);

eq3 = 23.5 == a*(1 - exp(b*59.5)) + c*(exp(d*59.5) - 1);

eq4 = 13.7 == a*(1 - exp(b*35.7)) + c*(exp(d*35.7) - 1);

[a,b,c,d]=vpasolve(eq1,eq2,eq3,eq4,a,b,c,d)

[John D'Errico]

That is what I would try myself. But I think vpasolve gets lost here.

[a,b,c,d]=vpasolve(eq1,eq2,eq3,eq4,a,b,c,d)
a =
475.95580519421968311347439860569
b =
-177.97484365049759450635900054123
c =
0
d =
28.872288438566491015788619039477

Since C is identically zero, the second term in each equation goes away.

a*(1-exp(b*[8.2 63 59.5 35.7]))
ans =
[ 475.95580519421968311347439860569, 475.95580519421968311347439860569, 475.95580519421968311347439860569, 475.95580519421968311347439860569]

The problem arises because, if you compute

b*[8.2 63 59.5 35.7]
ans =
[ -1459.3937179340802749521438044381, -11212.415149981348453900617034098, -10589.503197204606873128360532203, -6353.701918322764123877016319322]

Now, what happens when you raise e to each of those powers? The result is a massive underflow, that even vpasolve and the symbolic toolbox fails to get right.

[madhan ravi]

@ John D’Errico you are right thank you for pointing it out, how to achieve precision?

[John D'Errico]

By the way, +1 for your answer, because it is how one would normally solve the problem, and in fact, what I might have tried myself in a sleepy moment.

There is a LOT in your question. First, I would want to point out that this problem is essentially insolvable in double precision arithmetic.

For example, just look at the first two equations, and for now, completely ignore the second term. VPASOLVE did that any way. ;-)

eq1 = 3.8 == a*(1 - exp(b*8.2)) + c*(exp(d*8.1) - 1);
eq2 = 27.9 == a*(1 - exp(b*63)) + c*(exp(d*35.667) - 1);

Essentially, consider the terms

(1 - exp(b*8.2))
(1 - exp(b*63))

Now, give b any value. I'll pick -1 for b, just for a moment, since it is actually irrelevant for this discussion at this moment.

b = -1;
exp(b*8.2)
ans =
      0.000274653569972143

exp(b*63)
ans =
      4.35961000006308e-28

So you can see that exp(8.2*b) and exp(63*b) are different by roughly 25 powers of 10. Then, when we try to subtract those results from 1, we must get numerical garbage for at least one of them. The same thing happens in the second set of terms. Whatever the value of d, when compared to 1, exp(d*K), where K is a constant that varies between 8.1 and 59.5 is either a huge number, too may orders of magnitude larger than 1, or too many orders of magnitude smaller than 1.

So these equations are essentially garbage production facilities, when viewed in double precision arithmetic. No matter than the values of b and d are, they will contain uncomputable things in any way you try to use them. So while you MIGHT decide to try to solve them using syms, it won't matter because you could not actually verify the answer is correct in MATLAB. This is actually an important issue here. Why do I say that?

Suppose you could compute an answer that can be shown to be valid when using symbolic computations with hundreds of significant digits? In fact, the left hand side of those equations are known only to 2 significant digits. In effect, those numbers on the left hand side are rounded. You don't actually know the values those numbers really should have taken on. But the right hand sides, with those nasty exponentials will be highly ill-posed. Change the left hand side numbers by a tiny fraction, and the unknown coefficients {a,b,c,d} will potentially be impacted by a massive amount. And since we can't actually know the "true" value of those "rounded" numbers on the LHS of the equations, this problem becomes a truly difficult one, and probably one without a meaningful answer.

Realistically one might want to think about the problem in terms of interval arithmetic, just to see the difficulties involved. So, rather than the number 3.8, we should think about the interval [3.75,3.85]. Similarly, 8.2 and 8.1 are not those exact numbers, since they are just rounded forms. The problem now would look like (in terms of interval arithmetic, IF we had that capability in MATLAB) this:

eq1 = [3.75,3.85] == a*(1 - exp(b*[8.15,8.25])) + c*(exp(d*[8.05,8.15]) - 1);

Again, what I wrote is NOT valid MATLAB code to represent intervals for those numbers, but purely pseudo-code. The result, IF we could find one, would logically involve intervals for each of a,b,c,d. But, remember, that for any value of b, those exponentials blow up so large or get so small, that there will be serious numerical problems.

exp(1*[62.5 63.5])
ans =
      1.39124981295083e+27      3.78180908539129e+27

So whatever the parameters are, those exponentials will produce garbage in some form. The net result will be that the uncertainties in all of these rounded numbers will explode, to potentially make the solution have hugely wide uncertainties in the estimated parameters.

So, not only cannot you solve this problem using double precision, the solution is probably numerically meaningless anyway, since the supplied numbers are known to insufficient precision.

Ok, having said all that, suppose we knew that these numbers are all EXACTLY as they are given, that temporarily, we live in a universe where those numbers were all measured to have EXACTLY the values we are given. Could we then solve the problem? After all, in that case, a solution probably does exist, in THEORY, IF we were working in some hundreds of significant digits.

My first thought is don't even bother thinking you could replace those exponentials with truncated Taylor series. After all, if you could these equations would then reduce to approximate systems of polynomial equations. But remember that the sizes of those internal coefficients on b and d might require hundreds or more terms in the Taylor series to get convergence. So brute force will never suffice for something like this.

I'll end this comment here, to give me some time to think and play a bit in MATLAB....

[madhan ravi]

Highly appreciate your efforts to explain and it gave an idea to think logically and to suspect calculation in an effective way ?

[John D'Errico]

To reflect on the ill-posed nature of the problem, suppose we had the "true" values of b and d? I'll arbitrarily pick b=-1, d=-1 here.

digits 30
vpa(subs(eq1,[b,d],[-1 -1]))
vpa(subs(eq2,[b,d],[-1 -1]))
vpa(subs(eq3,[b,d],[-1 -1]))
vpa(subs(eq4,[b,d],[-1 -1]))

ans =
3.8 == 0.999725346430027857672372210606*a - 0.999696460861921133339134490468*c
ans =
27.9 == 0.999999999999999999999999999564*a - 0.999999999999999676392398964665*c
ans =
23.5 == 0.999999999999999999999999985563*a - 0.999999999999999999999999985563*c
ans =
13.7 == 0.99999999999999968689716782221*a - 0.99999999999999968689716782221*c

So we could write this in a partitioned iterative form, for any values of b and d, then try to solve for a and c. But see that the resulting linear system for a and c above will be an incredibly ill-conditioned system in itself. Pick any values of b and d, and those equations will change hugely, but always they will be massively ill-conditioned. That means the values of a and c will be incredibly sensitive to tiny perturbations in the values of b and d. Add in the fact that we really don't know the true values of those rounded constants, and we are left with a complete mess, something that would cause Cthulhu himself to run screaming from the room in terror. (An H.P. Lovecraft reference, in case anyone cares.)

[.999725346430027857672372210606,-0.999696460861921133339134490468
0.999999999999999999999999999564,-0.999999999999999676392398964665
0.999999999999999999999999985563,-0.999999999999999999999999985563
0.99999999999999968689716782221,-0.99999999999999968689716782221]

Again, rank or cond won't really do that matrix justice, if we try to do it in double precision, but I'll try.

cond([.999725346430027857672372210606,-0.999696460861921133339134490468
0.999999999999999999999999999564,-0.999999999999999676392398964665
0.999999999999999999999999985563,-0.999999999999999999999999985563
0.99999999999999968689716782221,-0.99999999999999968689716782221])
ans =
        159876.883630799

Hmm. So even cond suggests that for b=d=-1, any tiny perturbations in those unknown constants, all of which we know only to within about 2 significant digits, are then amplified by a factor of something on the order of roughly 160000 when we try to infer the values of a and c.

Other values of b and d will behave similarly nastily, but lets see what we can do, by following this path...

syms b d
M = [1 - exp(b*8.2), exp(d*8.1) - 1; ...
  1 - exp(b*63), exp(d*35.667) - 1; ...
  1 - exp(b*59.5), exp(d*59.5) - 1; ...
  1 - exp(b*35.7), exp(d*35.7) - 1];

B = [3.8;27.9;23.5;13.7];
ac = @(bd) pinv(vpa(subs(M,{b,d},bd),100))*B;

So, in fairly high precision, as a function of the unknown values of b and d, the function ac computes a and c.

ac([-1 -2])
ans =
 65194.8736654247620537283399937
 65173.1736654247522084970077505

ac([-1 -1])
ans =
 -619458.586882092496094139591758
 -619480.286882092539773390185907

ac([.5 0])
ans =
 -0.00000000000064823448769947172753851960683
                                            0

Again, what you need to recognize is just how sensitive a and c are to the values of b and d. I could probably use the above scheme to compute values for each of those parameters that fits the equations, for those specific constants. But they won't be of any real value, since we really don't know the values of those constants to a precision that is sufficient to estimate the parameters.

[madhan ravi]

Well explained @John D’Errico

[Walter Roberson]

My search routines have been unable to find a reasonable solution. I have identified an area with residue about 0.48, but it is pretty narrow. There are large areas where the residue is about 3.1 that are basically driving d to infinity to get rid of the c d terms, and if you do not happen to find the narrow better location then you can waste a lot of time on that other search region.

When I say narrow, I mean that I have found about 7 total points with residue less than 1. Not close ranges, individuals with numeric accidents.

Perhaps if I were to switch to symbolic

[madhan ravi]

After all your comments precision is becoming an essential thing to me now :) you people are really geniuses.

[Walter Roberson]

Looking again, there are a lot more points with very close residue than I had thought I had. Residues that are approximately 0.44011 were found for over 550 points, with a values ranging from 1337.8 to 171822.3, with b from -0.000388 to -2.965e-06, with c from 52.7 to 114.7, and d from -0.00244 to -0.00104 . a and b are practically linear with each other, and c and d are practically linear with each other, but the other pairs are only mostly linear with each other.

The residue was being calculated as

F = (a*(exp((41*b)/5) - 1) - c*(exp((81*d)/10) - 1) + 19/5)^2 + (a*(exp((119*b)/2) - 1) - c*(exp((119*d)/2) - 1) + 47/2)^2 + (a*(exp((357*b)/10) - 1) - c*(exp((357*d)/10) - 1) + 137/10)^2 + (a*(exp(63*b) - 1) - c*(exp((35667*d)/1000) - 1) + 279/10)^2

My tests did not show any sign of a "bowl" -- no minima that the values then rise again on the other side. I think doing better would require higher precision.

The marginally best combination I found was near 33762.2195025919 -1.50991408545558e-05 110.079764424126 -0.00109033033598055 but near 740.322188783004 -0.000712129376445178 37.1451062324705 -0.00368725773398294 is not far off

[madhan ravi]

Thank you once again @sir Walter.

Answer 2

[Alex Sha]

Here is the right solution, obtained from 1stOpt:

a: -2.38880243947768E-9
b: 0.357245116257183
c: -32.3471147348507
d: -0.0154282661155712

[Walter Roberson]

That does give a pretty good solution.

[Alex Sha]

fevl=

-5.85309578582383E-13

1.59872115546023E-12

-6.32383034826489E-12

-3.63442609341291E-12

want to know what is the “pretty good solution”?

[Walter Roberson]

The residue

F = (a*(exp((41*b)/5) - 1) - c*(exp((81*d)/10) - 1) + 19/5)^2 + (a*(exp((119*b)/2) - 1) - c*(exp((119*d)/2) - 1) + 47/2)^2 + (a*(exp((357*b)/10) - 1) - c*(exp((357*d)/10) - 1) + 137/10)^2 + (a*(exp(63*b) - 1) - c*(exp((35667*d)/1000) - 1) + 279/10)^2

is quite small, on the order of 1e-23 which is quite good for this situation.

[Alex Sha]

Hi, for "(a*(exp((41*b)/5) - 1) - c*(exp((81*d)/10) - 1) + 19/5)^2 + (a*(exp((119*b)/2) - 1) - c*(exp((119*d)/2) - 1) + 47/2)^2 + (a*(exp((357*b)/10) - 1) - c*(exp((357*d)/10) - 1) + 137/10)^2 + (a*(exp(63*b) - 1) - c*(exp((35667*d)/1000) - 1) + 279/10)^2", the minmum value：

Objective Function (Min.): 1.67130234831156E-23
a: -2.38880243946672E-9
b: 0.357245116257267
c: -32.347114734847
d: -0.015428266115563