Tuesday, 29 January 2013

Solving equations using SAGE

It is very easy to solve equations using the open source SAGE maths software. For example, to solve the equation:
we can type in SAGE:
> solve(3*x^2 + 16*x - 28 == 0, x)
[x == -2/3*sqrt(37) - 8/3, x == 2/3*sqrt(37) - 8/3]
We can round the answer to three decimal places by typing:
> round(-2/3*sqrt(37) - 8/3, 3)
-6.722
> round(2/3*sqrt(37) - 8/3, 3)    
1.389
This tells us that the two solutions are 1.389 and -6.722, rounded to three decimal places.
Easy peasy!

Another example
In the Little Book of R for Biomedical Statistics I have an equation:
n = numerator/denominator
where
numerator =  2 * ((qalpha + qgamma)^2) * pi0 * (1 - p0)
denominator = (piT - piC)^2
pi0 = (piT + piC)/2

Let's use SAGE to rearrange the equation in terms of piT: (Here I type 'a' for qalpha, 'g' for qgamma, 'T' for piT and 'C' for piC):
> a, g, T, C, n = var('a g T C n')
> numerator = 2 * ((a+g)^2) * ((T+C)/2) * (1 - ((T+C)/2))
> denominator = (T - C)^2
> solve(n == numerator/denominator, T)
[T == -((a^2 + 2*a*g + g^2 - 2*n)*C - a^2 - 2*a*g - g^2 + sqrt(-8*C^2*n + a^2 + 
2*a*g + g^2 + 8*C*n)*(a + g))/(a^2 + 2*a*g + g^2 + 2*n), T == -((a^2 + 2*a*g + 
g^2 - 2*n)*C - a^2 - 2*a*g - g^2 - sqrt(-8*C^2*n + a^2 + 2*a*g + g^2 + 
8*C*n)*(a + g))/(a^2 + 2*a*g + g^2 + 2*n)]

No comments: