MATH01, Linear equation systems
Back to the previous page
List of posts to read before reading this article
Contents
Equation
import sympy
#sympy.init_printing()
x = sympy.Symbol("x")
sympy.solve(x**2 + 2*x - 3)
OUTPUT
: \([-3, 1]\)
import sympy
#sympy.init_printing()
x = sympy.Symbol("x")
a, b, c = sympy.symbols("a, b, c")
sympy.solve(a * x**2 + b * x + c, x)
OUTPUT
: \(\left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]\)
import sympy
#sympy.init_printing()
x = sympy.Symbol("x")
sympy.solve(sympy.sin(x) - sympy.cos(x), x)
OUTPUT
: \(\left [ - \frac{3 \pi}{4}, \quad \frac{\pi}{4}\right ]\)
import sympy
#sympy.init_printing()
x = sympy.Symbol("x")
sympy.solve(sympy.exp(x) + 2 * x, x)
OUTPUT
: \(\left [ - \operatorname{LambertW}{\left (\frac{1}{2} \right )}\right ]\)
import sympy
#sympy.init_printing()
x = sympy.Symbol("x")
sympy.solve(x**5 - x**2 + 1, x)
OUTPUT
:
import sympy
sympy.init_printing()
x = sympy.Symbol("x")
y = sympy.Symbol("y")
eq1 = x + 2 * y - 1
eq2 = x - y + 1
sympy.solve([eq1, eq2], [x, y], dict=True)
OUTPUT
: \(\left [ \left \{ x : - \frac{1}{3}, \quad y : \frac{2}{3}\right \}\right ]\)
import sympy
sympy.init_printing()
x = sympy.Symbol("x")
y = sympy.Symbol("y")
eq1 = x**2 - y
eq2 = y**2 - x
sympy.solve([eq1, eq2], [x, y], dict=True)
SUPPLEMENT
import sympy
sympy.init_printing()
x = sympy.Symbol("x")
y = sympy.Symbol("y")
eq1 = x**2 - y
eq2 = y**2 - x
sols = sympy.solve([eq1, eq2], [x, y], dict=True)
[eq1.subs(sol).simplify() == 0 and eq2.subs(sol).simplify() == 0 for sol in sols]
OUTPUT
: [True, True, True, True]
OUTPUT
:
Square Systems
Main code : method1
from sympy import Matrix
#sympy.init_printing()
A = Matrix([[2,3],[5,4]])
b = Matrix([4,3])
x = A.solve(b)
print(x)
OUTPUT
: Matrix([[-1], [2]])
Main code : method1
import numpy as np
from scipy import linalg as la
A = np.array([[2,3],[5,4]])
b = np.array([4,3])
x = la.solve(A,b)
print(x)
OUTPUT
: [-1. 2.]
Rectangular Systems
Main code
from sympy import symbols, Matrix, solve
#sympy.init_printing()
x_vars = symbols("x_1, x_2, x_3")
A = Matrix([[1, 2, 3], [4, 5, 6]])
x = Matrix(x_vars)
b = Matrix([7, 8])
solution = solve(A*x - b, x_vars)
print(solution)
Data fitting
Main code
import numpy as np
from scipy import linalg as la
import matplotlib.pyplot as plt
# define true model parameters
x = np.linspace(-1, 1, 100)
a, b, c = 1, 2, 3
y_exact = a + b * x + c * x**2
# simulate noisy data
m = 100
X = 1 - 2 * np.random.rand(m)
Y = a + b * X + c * X**2 + np.random.randn(m)
# fit the data to the model using linear least square
A = np.vstack([X**0, X**1, X**2]) # see np.vander for alternative
sol, r, rank, sv = la.lstsq(A.T, Y)
y_fit = sol[0] + sol[1] * x + sol[2] * x**2
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(X, Y, 'go', alpha=0.5, label='Simulated data')
ax.plot(x, y_exact, 'k', lw=2, label='True value $y = 1 + 2x + 3x^2$')
ax.plot(x, y_fit, 'b', lw=2, label='Leat square fit')
ax.set_xlabel(r"$x$", fontsize=18)
ax.set_ylabel(r"$y$", fontsize=18)
ax.legend(loc=2)
plt.show()
Main code
import numpy as np
from scipy import linalg as la
import matplotlib.pyplot as plt
# define true model parameters
x = np.linspace(-1, 1, 100)
a, b, c = 1, 2, 3
y_exact = a + b * x + c * x**2
# simulate noisy data
m = 100
X = 1 - 2 * np.random.rand(m)
Y = a + b * X + c * X**2 + np.random.randn(m)
# fit the data to the model using linear least square:
# 1st order polynomial
A = np.vstack([X**n for n in range(2)])
sol, r, rank, sv = la.lstsq(A.T, Y)
y_fit1 = sum([s * x**n for n, s in enumerate(sol)])
# 15th order polynomial
A = np.vstack([X**n for n in range(16)])
sol, r, rank, sv = la.lstsq(A.T, Y)
y_fit15 = sum([s * x**n for n, s in enumerate(sol)])
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(X, Y, 'go', alpha=0.5, label='Simulated data')
ax.plot(x, y_exact, 'k', lw=2, label='True value $y = 1 + 2x + 3x^2$')
ax.plot(x, y_fit1, 'b', lw=2, label='Least square fit [1st order]')
ax.plot(x, y_fit15, 'm', lw=2, label='Least square fit [15th order]')
ax.set_xlabel(r"$x$", fontsize=18)
ax.set_ylabel(r"$y$", fontsize=18)
ax.legend(loc=2)
plt.show()
List of posts followed by this article
Reference