Solving System of Equations in Swift

Solving System of Equations in Swift

I do not know much about MATLAB, but A\b seems to return the solution of AX=B. For such purpose, you can use la_solve in LinearAlgebra framework (Accelerate.vecLib.LinearAlgebra) in iOS/macOS.

import Foundation
import Accelerate.vecLib.LinearAlgebra

let A: [Double] = [
1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 4.0, 1.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]
let matA = la_matrix_from_double_buffer(A, 10, 10, 10, la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))

let b: [Double] = [
0, -15, -15, -3, -3, 45, -12, -6, 0, 0
]
let vecB = la_matrix_from_double_buffer(b, 10, 1, 1, la_hint_t(LA_NO_HINT), la_attribute_t(LA_DEFAULT_ATTRIBUTES))

let vecCj = la_solve(matA, vecB)
var cj: [Double] = Array(repeating: 0.0, count: 10)

let status = la_matrix_to_double_buffer(&cj, 1, vecCj)
if status == la_status_t(LA_SUCCESS) {
print(cj) //->[0.0, -2.9185349611542728, -3.3258601553829079, 1.2219755826859044, -4.5620421753607099, 14.026193118756936, -6.5427302996670358, 0.14472807991120964, -0.036182019977802411, 0.0]
} else {
print("Failure: \(status)")
}

The result seems to be the same (other than the precision) as yours.

Solve Equations Of Type A*X = B Using dgtsv_ or sgtsv_

dgtsv_() expects the lower/main/upper diagonal of the tri-diagonal
matrix as separate arguments. You can pass the address of variable arrays with &.

All integer parameters are addresses of __CLPK_integer aka Int32
variables.

The right-hand side vector b is overwritten with the solution x to
the A x = b equation. The three vectors describing A are overwritten
as well, so you might want to make copies of the original data.

Example:

import Swift
import Accelerate

var mainDiagA = [ 1.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 1.0 ]
var upperDiagA = [ 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ]
var lowerDiagA = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0 ]

var b = [0.0, -15.0, -15.0, -3.0, -3.0, 45.0, -12.0, -6.0, 0.0, 0.0 ]

var n = Int32(mainDiagA.count) // Order of matrix A
var nrhs = Int32(1) // Number of right-hand sides
var info = Int32(0) // Result code

dgtsv_(&n, &nrhs, &lowerDiagA, &mainDiagA, &upperDiagA, &b, &n, &info)
if info == 0 { // success
print(b)
// [0.0, -2.9185349611542732, -3.3258601553829075, 1.2219755826859044, -4.5620421753607099, 14.026193118756938, -6.5427302996670367, 0.14472807991120964, -0.036182019977802411, 0.0]

}

Swift: Multiplication and brackets calculation doesn't work

This should work:

let numerator: Double = (5-2) * (10-5)
let denumerator: Double = (4-2) * (10-5)

Fist you calculate the numerator and denumerator. And finally the result:

print(result)
let result: Double = numerator/denumerator
//1.5

Generating a simple algebraic expression in swift

Since all you want is a string representing an equation and a value for x, you don't need to do any solving. Just start with x and transform it until you have a nice equation. Here's a sample: (copy and paste it into a Playground to try it out)

import UIKit

enum Operation: String {
case addition = "+"
case subtraction = "-"
case multiplication = "*"
case division = "/"

static func all() -> [Operation] {
return [.addition, .subtraction, .multiplication, .division]
}

static func random() -> Operation {
let all = Operation.all()
let selection = Int(arc4random_uniform(UInt32(all.count)))
return all[selection]
}

}

func addNewTerm(formula: String, result: Int) -> (formula: String, result: Int) {
// choose a random number and operation
let operation = Operation.random()
let number = chooseRandomNumberFor(operation: operation, on: result)
// apply to the left side
let newFormula = applyTermTo(formula: formula, number: number, operation: operation)
// apply to the right side
let newResult = applyTermTo(result: result, number: number, operation: operation)
return (newFormula, newResult)
}

func applyTermTo(formula: String, number:Int, operation:Operation) -> String {
return "\(formula) \(operation.rawValue) \(number)"
}

func applyTermTo(result: Int, number:Int, operation:Operation) -> Int {
switch(operation) {
case .addition: return result + number
case .subtraction: return result - number
case .multiplication: return result * number
case .division: return result / number
}
}

func chooseRandomNumberFor(operation: Operation, on number: Int) -> Int {
switch(operation) {
case .addition, .subtraction, .multiplication:
return Int(arc4random_uniform(10) + 1)
case .division:
// add code here to find integer factors
return 1
}
}

func generateFormula(_ numTerms:Int = 1) -> (String, Int) {
let x = Int(arc4random_uniform(10))
var leftSide = "x"
var result = x

for i in 1...numTerms {
(leftSide, result) = addNewTerm(formula: leftSide, result: result)
if i < numTerms {
leftSide = "(" + leftSide + ")"
}
}

let formula = "\(leftSide) = \(result)"

return (formula, x)
}

func printFormula(_ numTerms:Int = 1) {
let (formula, x) = generateFormula(numTerms)
print(formula, " x = ", x)
}

for i in 1...30 {
printFormula(Int(arc4random_uniform(3)) + 1)
}

There are some things missing. The sqrt() function will have to be implemented separately. And for division to be useful, you'll have to add in a system to find factors (since you presumably want the results to be integers). Depending on what sort of output you want, there's a lot more work to do, but this should get you started.

Here's sample output:

(x + 10) - 5 = 11                       x =  6
((x + 6) + 6) - 1 = 20 x = 9
x - 2 = 5 x = 7
((x + 3) * 5) - 6 = 39 x = 6
(x / 1) + 6 = 11 x = 5
(x * 6) * 3 = 54 x = 3
x * 9 = 54 x = 6
((x / 1) - 6) + 8 = 11 x = 9

Swift - Quadratic equation solver gives wrong results?

You erroneously multiply by a instead of dividing

( -b + sqrt(delta)) / (2*a)

Function which returns the least-squares solution to a linear matrix equation

The Accelerate framework included the LAPACK linear algebra package,
which has a DGELS function to solve under- or overdetermined linear systems. From the documentation:

DGELS solves overdetermined or underdetermined real linear systems
involving an M-by-N matrix A, or its transpose, using a QR or LQ
factorization of A. It is assumed that A has full rank.

Here is an example how that function can be used from Swift.
It is essentially a translation of this C sample code.

func solveLeastSquare(A A: [[Double]], B: [Double]) -> [Double]? {
precondition(A.count == B.count, "Non-matching dimensions")

var mode = Int8(bitPattern: UInt8(ascii: "N")) // "Normal" mode
var nrows = CInt(A.count)
var ncols = CInt(A[0].count)
var nrhs = CInt(1)
var ldb = max(nrows, ncols)

// Flattened columns of matrix A
var localA = (0 ..< nrows * ncols).map {
A[Int($0 % nrows)][Int($0 / nrows)]
}

// Vector B, expanded by zeros if ncols > nrows
var localB = B
if ldb > nrows {
localB.appendContentsOf([Double](count: ldb - nrows, repeatedValue: 0.0))
}

var wkopt = 0.0
var lwork: CInt = -1
var info: CInt = 0

// First call to determine optimal workspace size
dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows, &localB, &ldb, &wkopt, &lwork, &info)
lwork = Int32(wkopt)

// Allocate workspace and do actual calculation
var work = [Double](count: Int(lwork), repeatedValue: 0.0)
dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows, &localB, &ldb, &work, &lwork, &info)

if info != 0 {
print("A does not have full rank; the least squares solution could not be computed.")
return nil
}
return Array(localB.prefix(Int(ncols)))
}

Some notes:

  • dgels_() modifies the passed matrix and vector data, and expects
    the matrix as "flat" array containing the columns of A.
    Also the right-hand side is expected as an array with length max(M, N).
    For this reason, the input data is copied to local variables first.
  • All arguments must be passed by reference to dgels_(), that's why
    they are all stored in vars.
  • A C integer is a 32-bit integer, which makes some conversions between
    Int and CInt necessary.

Example 1: Overdetermined system, from http://www.seas.ucla.edu/~vandenbe/103/lectures/ls.pdf.

let A = [[ 2.0, 0.0 ],
[ -1.0, 1.0 ],
[ 0.0, 2.0 ]]
let B = [ 1.0, 0.0, -1.0 ]
if let x = solveLeastSquare(A: A, B: B) {
print(x) // [0.33333333333333326, -0.33333333333333343]
}

Example 2: Underdetermined system, minimum norm
solution to x_1 + x_2 + x_3 = 1.0.

let A = [[ 1.0, 1.0, 1.0 ]]
let B = [ 1.0 ]
if let x = solveLeastSquare(A: A, B: B) {
print(x) // [0.33333333333333337, 0.33333333333333337, 0.33333333333333337]
}

Update for Swift 3 and Swift 4:

func solveLeastSquare(A: [[Double]], B: [Double]) -> [Double]? {
precondition(A.count == B.count, "Non-matching dimensions")

var mode = Int8(bitPattern: UInt8(ascii: "N")) // "Normal" mode
var nrows = CInt(A.count)
var ncols = CInt(A[0].count)
var nrhs = CInt(1)
var ldb = max(nrows, ncols)

// Flattened columns of matrix A
var localA = (0 ..< nrows * ncols).map { (i) -> Double in
A[Int(i % nrows)][Int(i / nrows)]
}

// Vector B, expanded by zeros if ncols > nrows
var localB = B
if ldb > nrows {
localB.append(contentsOf: [Double](repeating: 0.0, count: Int(ldb - nrows)))
}

var wkopt = 0.0
var lwork: CInt = -1
var info: CInt = 0

// First call to determine optimal workspace size
var nrows_copy = nrows // Workaround for SE-0176
dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows_copy, &localB, &ldb, &wkopt, &lwork, &info)
lwork = Int32(wkopt)

// Allocate workspace and do actual calculation
var work = [Double](repeating: 0.0, count: Int(lwork))
dgels_(&mode, &nrows, &ncols, &nrhs, &localA, &nrows_copy, &localB, &ldb, &work, &lwork, &info)

if info != 0 {
print("A does not have full rank; the least squares solution could not be computed.")
return nil
}
return Array(localB.prefix(Int(ncols)))
}

Generating random doable math problems swift

I started programming with Turbo Pascal and as Niklaus Wirth said: Algorithms + Data Structure = Programs. You need to define data structures appropriate for your program.

First, some basic data structures. (Swift enum is much more powerful than that in other languages)

enum MathElement : CustomStringConvertible {
case Integer(value: Int)
case Percentage(value: Int)
case Expression(expression: MathExpression)

var description: String {
switch self {
case .Integer(let value): return "\(value)"
case .Percentage(let percentage): return "\(percentage)%"
case .Expression(let expr): return expr.description
}
}

var nsExpressionFormatString : String {
switch self {
case .Integer(let value): return "\(value).0"
case .Percentage(let percentage): return "\(Double(percentage) / 100)"
case .Expression(let expr): return "(\(expr.description))"
}
}
}

enum MathOperator : String {
case plus = "+"
case minus = "-"
case multiply = "*"
case divide = "/"

static func random() -> MathOperator {
let allMathOperators: [MathOperator] = [.plus, .minus, .multiply, .divide]
let index = Int(arc4random_uniform(UInt32(allMathOperators.count)))

return allMathOperators[index]
}
}

Now the main class:

class MathExpression : CustomStringConvertible {
var lhs: MathElement
var rhs: MathElement
var `operator`: MathOperator

init(lhs: MathElement, rhs: MathElement, operator: MathOperator) {
self.lhs = lhs
self.rhs = rhs
self.operator = `operator`
}

var description: String {
var leftString = ""
var rightString = ""

if case .Expression(_) = lhs {
leftString = "(\(lhs))"
} else {
leftString = lhs.description
}
if case .Expression(_) = rhs {
rightString = "(\(rhs))"
} else {
rightString = rhs.description
}

return "\(leftString) \(self.operator.rawValue) \(rightString)"
}

var result : Any? {
let format = "\(lhs.nsExpressionFormatString) \(`operator`.rawValue) \(rhs.nsExpressionFormatString)"
let expr = NSExpression(format: format)
return expr.expressionValue(with: nil, context: nil)
}

static func random() -> MathExpression {
let lhs = MathElement.Integer(value: Int(arc4random_uniform(10)))
let rhs = MathElement.Integer(value: Int(arc4random_uniform(10)))

return MathExpression(lhs: lhs, rhs: rhs, operator: .random())
}
}

Usage:

let a = MathExpression(lhs: .Integer(value: 1), rhs: .Integer(value: 2), operator: .divide)
let b = MathExpression(lhs: .Integer(value: 10), rhs: .Percentage(value: 20), operator: .minus)
let c = MathExpression.random()

let d = MathExpression(lhs: .Integer(value: 1), rhs: .Integer(value: 2), operator: .plus)
let e = MathExpression(lhs: .Integer(value: 3), rhs: .Integer(value: 4), operator: .plus)
let f = MathExpression(lhs: .Expression(expression: d), rhs: .Expression(expression: e), operator: .multiply)

let x = MathExpression.random()
let y = MathExpression.random()
let z = MathExpression(lhs: .Expression(expression: x), rhs: .Expression(expression: y), operator: .plus)

print("a: \(a) = \(a.result!)")
print("b: \(b) = \(b.result!)")
print("c: \(c) = \(c.result!)")

print("f: \(f) = \(f.result!)") // the classic (1 + 2) * (3 + 4)
print("z: \(z) = \(z.result!)") // this is completely random

csparse error while solving sparse system of equations

Matrix A above is in matrix tripplet form. First we need to convert it to compressed column format (B) and then apply the cs_chsol function to get the results.

// Convert matrix tripplet to column compressed form

cs *B = cs_compress(&A);
int status = cs_cholsol(0,&B,&b[0]);

I implemented it in my program and everything is running perfectly fine now.

Numerical Solution of System of Difference/Differential Algebraic Equations in Maxima

Well, if the goal is to compute a solution for the equations given, I think you can do that with so-called memoizing functions (i.e., functions which compute a result and remember it). In Maxima such functions are defined by f[k] := ... (instead of f(k) := ... for ordinary, non-memoizing functions). In this case I think you would have:

Y[t] := A[t]*K[t]^(1/3)*Ly[t]^(2/3);
K[t] := K[t - 1] + s*Y[t - 1] - d*K[t - 1];
A[t] := A[t - 1] + z*A[t - 1]*La[t - 1];
Ly[t] := L - La[t];
La[t] := l*L;

Y[0] : 9663.8253;
K[0] : 100;
A[0] : 100;
Ly[0] : 95.0;
La[0] : 5.0;

s: 0.15;
d: 0.07;
z: 0.02;
l: 0.05;
L: 100.0;

and then you can compute e.g. makelist ([Y[n], K[n], A[n], Ly[n], La[t]], n, 1, 10);

That works if you can rearrange the equations yourself; in this case it's easy. If you need to solve other problems and that approach fails, you'll need something more powerful. I don't know what to suggest. My advice is to ask on the maxima-discuss@lists.sourceforge.net mailing list where most discussion about Maxima takes place.



Related Topics



Leave a reply



Submit