how to solve these cubic equations with three terms and three equations

syms p a b
E1 = 56.*p + (1-p).*56.*exp(-a.*1.03699)-56.*exp(-b.*1.03699) == 0;
E2 = 57.*p + (1-p).*57.*exp(-a.*1.05554)-57.*exp(-b.*1.05554) == 0;
E3 = 58.*p + (1-p).*58.*exp(-a.*1.07404)-58.*exp(-b.*1.07404) == 0;
result = solve(E1,E2,E3);

I want to know values of a, p , b

but the results indicate that a,p,b are 0, What should I do?

Answer 1

[Walter Roberson]

Q = @(v) sym(v);
syms p a b
E1 = 56.*p + (1-p).*56.*exp(-a.*Q(1.03699))-56.*exp(-b.*Q(1.03699)) == 0;
E2 = 57.*p + (1-p).*57.*exp(-a.*Q(1.05554))-57.*exp(-b.*Q(1.05554)) == 0;
E3 = 58.*p + (1-p).*58.*exp(-a.*Q(1.07404))-58.*exp(-b.*Q(1.07404)) == 0;
sol = vpasolve([E1, E2, E3], randn(3,1))

sol = struct with fields:
    a: 0.54360320751138251896224670429114
    b: 1.8957551477198055480839774519014e-39
    p: 1.0

sol = vpasolve([E1, E2, E3], randn(3,1))

sol = struct with fields:
    a: -1.5310493253716739735043666593253
    b: -2.7320581392130059765958037120213e-39
    p: 1.0

sol = vpasolve([E1, E2, E3], randn(3,1))

sol = struct with fields:
    a: -2.1128525577424056479951278913567
    b: -1.8655458644816284351600623851951e-39
    p: 1.0

sol = vpasolve([E1, E2, E3], randn(3,1))

sol = struct with fields:
    a: -0.41938339355777769244433988252232
    b: -0.41938339355777769244433988252232
    p: 8.4696183680066730034095698192125e-38

... when p = 1 then 1-p is 0 and the a term vanishes. This naturally leads to b = 0.

[John D'Errico]

@Walter Roberson

The solution found by you is that effectively, only a trivial solution seems to exist.

One question might be if this is an artifact of the use of floating point numbers, of the form 1.03699, so only 6 decimal digit approximations. At first I wondered if that may be the issue.

But then, looking carefully at the equations, we see that each of them are of the same identical form

E1 = 56.*p + (1-p).*56.*exp(-a.*1.03699)-56.*exp(-b.*1.03699) == 0;

Note that the right hand side is 0. So we can factor out (and divide by) the constant coefficient 56 from the first equation. Similarly, we can factor out, and divide by 57 for the second equation, and then also divde the third equation by 58. That leaves us with 3 somewhat simpler equations.

Q = @(v) sym(v);
syms p a b
E1 = p + (1-p).*exp(-a.*Q(1.03699))-exp(-b.*Q(1.03699)) == 0;
E2 = p + (1-p).*exp(-a.*Q(1.05554))-exp(-b.*Q(1.05554)) == 0;
E3 = p + (1-p).*exp(-a.*Q(1.07404))-exp(-b.*Q(1.07404)) == 0;

My guess is you are correct, in that there is only one (trivial) solution. Those three equations look too similar to each other that I would postulate there is no non-trivial solution.

[재혁]

Thanks for answer my question. but I want to like to obtain three distinct roots. Even though it might be feasible manually by taking the logarithm of both sides. If I use your method, I'm also curious about how to set the initial guess x0. Thanks.

Answer 2

[Sam Chak]

Annyeonghaseyo @재혁

Can you verify if the following three simplified equations are 100% correct, or are the equalities merely approximations?

Have you ever heard of half-life (not the popular first-person shooter game)? In nuclear physics, half-life refers to the amount of time it takes for half of the radioactive atoms in a sample to undergo decay.

By definition, radioactive decay is given by the equation:

where is the half-life of the radioactive element.

However, the equation can also be rearranged as:

where the symbol represents the decay constant (not associated with the logo of the first-person shooter game).

Notice the exponential function? Translating that into mathematical language, it means that the total mass of radioactive isotopes never actually reaches zero in finite time.

Back to your original question, the problem can be formulated to solve this equation:

so that .

Before we solve the problem, we can actually guess that either or . If we assume the solution , then we need to solve:

, which can be rearranged to:

, and it reduces to:

This means that the symbol a is a free parameter, and b equals it.

The non-zero solution for p has been demonstrated by @Walter Roberson and explained by @John D'Errico.

[Sam Chak]

@재혁, Let's test this out: 😅

p   = 0;
a   = 1;    % <-- freely assign a real value to this parameter
b   = a;
Eq1 = 56*p + (1 - p)*56*exp(- a*1.03699) - 56*exp(- b*1.03699)

Eq1 = 0

Eq2 = 57*p + (1 - p)*57*exp(- a*1.05554) - 57*exp(- b*1.05554)

Eq2 = 0

Eq3 = 58*p + (1 - p)*58*exp(- a*1.07404) - 58*exp(- b*1.07404)

Eq3 = 0