A very nonlinear system of three equations

asked 2013-09-05 04:15:08 +0200

OrionNav
1 ●1 ●1 ●2

updated 2023-05-19 22:13:40 +0200

5125 ●3 ●43 ●111

Here's a fun little problem: determine the exponential curve f(x) = c + ba^x defined by three points, (2,10), (4,6), and (5,5).

The system of three equations and three unknowns is

10 = c + ba^2

6 = c + ba^4

5 = c + ba^5

It's not that hard to solve numerically. With a little algebraic substitution and iteration, the answer turns out to be

a = 0.640388203

b = 16.53456516

c = 3.219223594

But is there a more elegant way to use Sage to arrive at this result? I'm stumped.

edit retag flag offensive close merge delete

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2013-09-05 16:42:10 +0200

Mark
106 ●2 ●5 ●11

If you just want to have a numerical solution to a small set of (nonlinear) equations, why not simply 'fsolve' them?

from scipy.optimize import fsolve

def equations(p):
    a, b, c = p
    return (c + b*abs(a)**2 - 10, c + b*abs(a)**4 - 6, c + b*abs(a)**5 - 5)

a,b,c = fsolve(equations, (1, 1, 1))

print a,b,c
print equations((a,b,c))
0.640388203202 16.5345651552 3.2192235936
(-2.5757174171303632e-13, -2.3625545964023331e-13,-1.9895196601282805e-13)

The 'abs()' is there to provide at least some additional constraining on a,b,c

I don't know if this is elegant - but it is straightforward.

edit flag offensive delete link

Comments

Thanks. I suspect that some or all of this is Python code, and I know absolutely nothing about Python. (I'm a decent Perl coder, though.) The Sage Reference Manual has information about maple.fsolve and qfsolve, but nothing about plain old fsolve.

OrionNav ( 2013-09-06 00:34:02 +0200 )edit

add a comment

answered 2013-09-05 09:23:45 +0200

kcrisman
12222 ●41 ●135 ●254

Have you tried find_fit? Presumably with an exact fit this would work very nicely. (Let me know if you need more details or if that doesn't work.)

That said, this isn't "elegant".

edit flag offensive delete link

Comments

Thanks! That works well enough. My first attempt gave an unexpected answer: model(x) = c+b*a^x find_fit([(2,10), (4,6), (5,5)], model) [a == -0.39038820320220763, b == 30.965434844830888, c == 5.2807764064044145] which is correct as far as it goes, but raising -0.39 to non-integer values of x results in a complex number. Giving it a non-default initial guess sussed out the expected answer: model(x) = c+b*a^x find_fit([(2,10), (4,6), (5,5)], model, initial_guess=(1,1,3)) [a == 0.6403882032022076, b == 16.534565155169098, c == 3.219223593595585] and I'm pleased with the accuracy of the result. Since this can be done with only two lines of code, I consider it more elegant than the fsolve approach (above).

OrionNav ( 2013-09-06 00:27:39 +0200 )edit

Yes, sometimes you do need to give it a little help! That's typical of such numerical routines.

kcrisman ( 2013-09-06 11:49:23 +0200 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

A very nonlinear system of three equations

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

A very nonlinear system of three equations edit

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

A very nonlinear system of three equations