CODEWITHSUNDEEP
pythonalgorithmmatlabnumpymath
plywoods
plywoods

Reputation: 257

How to resolve function approximation task in Python?

Consider the complex mathematical function on the line [1, 15]: f(x) = sin(x / 5) * exp(x / 10) + 5 * exp(-x / 2)

enter image description here

polynomial of degree n (w_0 + w_1 x + w_2 x^2 + ... + w_n x^n) is uniquely defined by any n + 1 different points through which it passes. This means that its coefficients w_0, ... w_n can be determined from the following system of linear equations:

enter image description here

Where x_1, ..., x_n, x_ {n + 1} are the points through which the polynomial passes, and by f (x_1), ..., f (x_n), f (x_ {n + 1}) - values that it must take at these points.

I'm trying to form a system of linear equations (that is, specify the coefficient matrix A and the free vector b) for the polynomial of the third degree, which must coincide with the function f at points 1, 4, 10, and 15. Solve this system using the scipy.linalg.solve function.

A = numpy.array([[1., 1., 1., 1.], [1., 4., 8., 64.], [1., 10., 100., 1000.], [1., 15., 225., 3375.]])

V = numpy.array([3.25, 1.74, 2.50, 0.63])

numpy.linalg.solve(A, V)

I got the wrong answer, which isenter image description here

So the question is: is the matrix correct?

Upvotes: 3

Views: 9009

Answers (2)

Vadim Kirilchuk
Vadim Kirilchuk

Reputation: 3542

Is it from data science course? :) Here is an almost generic solution I did:

%matplotlib inline

import numpy as np;
import math;
import matplotlib.pyplot as plt;

def f(x):
    return np.sin(x / 5) * np.exp(x / 10) + 5 * np.exp(-x / 2)

# approximate at the given points (feel free to experiment: change/add/remove)
points = np.array([1, 4, 10, 15])
n = points.size

# fill A-matrix, each row is 1 or xi^0, xi^1, xi^2, xi^3 .. xi^n
A = np.zeros((n, n))
for index in range(0, n):
    A[index] = np.power(np.full(n, points[index]), np.arange(0, n, 1))

# fill b-matrix, i.e. function value at the given points
b = f(points)

# solve to get approximation polynomial coefficents
solve = np.linalg.solve(A,b)

# define the polynome approximation of the function
def polinom(x): 
    # Yi = solve * Xi where Xi = x^i
    tiles = np.tile(x, (n, 1))
    tiles[0] = np.ones(x.size)
    for index in range(1, n):
        tiles[index] = tiles[index]**index
    return solve.dot(tiles)

# plot the graphs of original function and its approximation
x = np.linspace(1, 15, 100)
plt.plot(x, f(x))
plt.plot(x, polinom(x))

# print out the coefficients of polynome approximating our function
print(solve)

Upvotes: 0

Rory Daulton
Rory Daulton

Reputation: 22544

No, your matrix is not correct.

The biggest mistake is your second sub-matrix for A. The third entry should be 4**2 which is 16 but you have 8. Less important, you have only two decimal places for your constants array V but you really should have more precision than that. Systems of linear equations are sometimes very sensitive to the provided values, so make them as precise as possible. Also, the rounding in your final three entries is bad: you rounded down but you should have rounded up. If you really want two decimal places (which I do not recommend) the values should be

V = numpy.array([3.25, 1.75, 2.51, 0.64])

But better would be

V = numpy.array([3.252216865271419, 1.7468459495903677,
                 2.5054164070002463, 0.6352214195786656])

With those changes to A and V I get the result

array([ 4.36264154, -1.29552587,  0.19333685, -0.00823565])

I get these two sympy plots, the first showing your original function and the second using the approximated cubic polynomial.

enter image description here

enter image description here

They look close to me! When I calculate the function values at 1, 4, 10, and 15, the largest absolute error is for 15, namely -4.57042132584462e-6. That is somewhat larger than I would have expected but probably is good enough.

Upvotes: 3

Related Questions