Solving a simultaneous equation through code

Question

This seems like an incredibly simple and silly question to ask, but everything I've found about it has been too complex for me to understand.

I have these two very basic simultaneous equations:

X = 2x + 2z
Y = z - x

Given that I know both X and Y, how would I go about finding x and z? It's very easy to do it by hand, but I have no idea how this would be done in code.

Answer 1

Other answers are missing the non-unique solution cases, and will silently fail on those cases. Since they're perfectly valid cases, it will likely help to check first.

Here's a heavily commented example. It's in Gleam, since the declarative style and explicit types help with illustration, and I hope is self-explanatory for non-Gleam speakers. To aid in understanding, note:

• /. is just the floating point version of /
• { and } are used where parentheses are usually used for order of operations reasons
• the case statement checks both parallel and coincident, and _ means "don't care".
• #(x, y) is a tuple with elements x and y.

The function proceeds by determining first whether the two equations represent lines that are parallel or not, and if so, coincident or not. It then splits into the three cases. If the lines are not parallel, the solution is unique. If they are parallel, there are either infinite solutions all the way along the line if the lines are coincident, or zero solutions if they are not.

For the unique solution case it uses Cramer's rule, which is derived in the comments using the elimination method. See other answers for why this may be undesirable.

pub type TwoByTwoSolution {
  Unique(#(Float,Float))   // 1 solution:  #(x, y) where the sol is x=x, y=y
  Infinite(#(Float,Float)) // ∞ solutions: #(a, b) where sols are x=t, y=a+bt for any t
  None                     // 0 solutions
}

// Solve two simultaneous equations in two unknowns, where:
//   a1*x + b1*y = c1 (1)
//   a2*x + b2*y = c2 (2)
// Zero-value coefficients are not supported!
pub fn solve_two_by_two_equations(a1: Float, b1: Float, c1: Float,
                                  a2: Float, b2: Float, c2: Float) -> TwoByTwoSolution
{
  //Possible floating point imprecision here. To evaluate impact,
  //consider the floating point precision of the x calculation below.
  let parallel   = a1/.a2 == b1/.b2
  let coincident = parallel && b1/.b2 == c1/.c2

  case parallel, coincident {

    False, _      -> {
      // (1) * b2  => a1*x*b2 + b1*y*b2 = c1*b2  ... (3)
      // (2) * b1  => a2*x*b1 + b2*y*b1 = c2*b1  ... (4)
      // (3) - (4) => a1*x*b2 - a2*x*b1 = c1*b2 - c2*b1
      //           => x*(a1*b2 - a2*b1) = c1*b2 - c2*b1
      //           => x = (c1*b2 - c2*b1) / (a1*b2 - a2*b1)
      let x = {c1*.b2 -. c2*.b1} /. {a1*.b2 -. a2*.b1}
      // (1)       => b1*y = c1 - a1*x
      //           =>    y = (c1 - a1*x) / b1
      Unique(#(x, {c1 -. a1*.x} /. b1))
    }

    True,  False  -> {
      None
    }

    True,  True   -> {
      // (1),x=t  => y = (c1 - a1*t) / b1
      //           = c1/b1 + (-a1/b1)*t
      Infinite(#(c1/.b1, -1.*.a1/.b1))
    }

  }
}

Answer 2

This seems like an incredibly simple and silly question to ask

Not at all. This is a very good question, and it has unfortunately a complex answer. Let's solve

a * x + b * y = u
c * x + d * y = v

I stick to the 2x2 case here. More complex cases will require you to use a library.

The first thing to note is that Cramer formulas are not good to use. When you compute the determinant

a * d - b * c

as soon as you have a * d ~ b * c, then you have catastrophic cancellation. This case is typical, and you must guard against it.

The best tradeoff between simplicity / stability is partial pivoting. Suppose that |a| > |c|. Then the system is equivalent to

a * c/a * x + bc/a * y = uc/a
      c * x +    d * y = v 

which is

cx + bc/a * y = uc/a
cx +       dy = v  

and now, substracting the first to the second yields

cx +       bc/a * y = uc/a
     (d - bc/a) * y = v - uc/a

which is now straightforward to solve: y = (v - uc/a) / (d - bc/a) and x = (uc/a - bc/a * y) / c. Computing d - bc/a is stabler than ad - bc, because we divide by the biggest number (it is not very obvious, but it holds -- do the computation with very close coefficients, you'll see why it works).

Now, if |c| > |a|, you just swap the rows and proceed similarly.

In code (please check the Python syntax):

def solve(a, b, c, d, u, v):
    if abs(a) > abs(c):
         f = u * c / a
         g = b * c / a
         y = (v - f) / (d - g)
         return ((f - g * y) / c, y)
    else:
         f = v * a / c
         g = d * a / c
         x = (u - f) / (b - g)
         return (x, (f - g * x) / a)

You can use full pivoting (requires you to swap x and y so that the first division is always by the largest coefficient), but this is more cumbersome to write, and almost never required for the 2x2 case.

For the n x n case, all the pivoting stuff is encapsulated into the LU decomposition, and you should use a library for this.

Answer 3

Algorithm to sove: Ax + By = C, Dx + Ey = F

START

1)  PRINT "Enter A"
2)  INPUT A
3)  PRINT "Enter B"
4)  INPUT B
5)  PRINT "Enter C"
6)  INPUT C
7)  PRINT "Enter D"
8)  INPUT D
9)  PRINT "Enter E"
10) INPUT E
11) PRINT "Enter F"
21) INPUT F

22) detS <-- AE - BC
22) detX <-- CE - BF
24) detY <-- AF - CD

25) x <-- detX / detS
26) y <-- detY / detS

27) PRINT x, y

END

Answer 4

@Alexandre , you missed one condition.. Here is the final code

void SolveLinearEquations (float a,float b,float c,float d,float u,float v, float &x, float &y)
{
float f;
float g;
if (abs(a) > abs(c))
{
     f = u * c / a;
     g = b * c / a;
     y = (v - f) / (d - g);
     if(c != 0)
        x = (f - g * y) / c;
     else
        x = (u - b * y)/a;
}
else
{
    f = v * a / c;
     g = d * a / c;
     x = (u - f) / (b - g);
     if (a != 0)
        y = (f - g * x) / a ;
     else
        y = (v - d * x)/c;
}
}

Answer 5

The following function might be useful for some:

function solve(s1,s2){    //only works if coefficients > 0

   str=s1 + " " +s2
   str=str.replace(/[^0123456789.-]/g, ' ') //eliminate letters
   str=str.replace( /\s\s+/g, ' ' )         //no double spaces


   var n=str.split(" ");            //put into an array
   var a=0,b=1,c=2,d=3,e=4,f=5      //see what were doing
   var x = ( n[c]*n[e] -n[b]*n[f])/(n[a]*n[e] - n[b]*n[d])
   var y= (n[c]-n[a]*x)/n[b]

   return({x:x, y:y})
}

To use:

result=solve("12x +  2y =32", "9x -5y=55")
alert (result.x+" ----- "+result.y)

Answer 6

(1)  ax + by = c    
(2)  dx + dy = f

(3)1*d     adx + bdy  = cd
(4)2*b     abx + bdy  = fb

(3)-(4)   adx - abx  = cd - fb

          x(ad-ab) = cd - fb


          x = (c*d - f*b)/(a*d-a*b) //use this equation for x



ax + by = c
     by = c - ax

     y = (c - a*x)/b    //use this equation for y