-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathLagrange.py
More file actions
66 lines (50 loc) · 1.96 KB
/
Lagrange.py
File metadata and controls
66 lines (50 loc) · 1.96 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 30 18:11:49 2021
@author: moatasim
"""
from sympy import *
from sympy import init_printing
import numpy as np
import matplotlib.pyplot as plt
print("This script constructs the Lagrange Polynomial from data generated by a function in a specified range. Note: when inputting the function, use ** to show presence of exponent, followed by the number or variable in the power. Also, input capital E for euler's number instead of just e")
print("WARNING: when inputting function, do not use quotation marks, otherwise sympify is not able to parse the input string!")
x=symbols('x')
print("for quick checking and to avoid errors, copy paste this= (E**x)+(2**(-x))+2*cos(x) - 6")
inputfunc = input("Enter the function: ")
func=sympify(inputfunc)
print("Since the question asks only for data points in [-2,2], we will not be asking user to input start and end points, we will just simply take [-2,2]")
datapoints=np.linspace(-2,2,10,endpoint=False)
print("Our set of x values: ", datapoints)
fx=[]
def calcfunc(a):
return func.evalf(subs={x:a})
for val in datapoints:
fx.append(calcfunc(val))
print("Our set of y values: ", fx)
print("By now we have constructed two lists, containing x and y values (our DATA), and the lagrange polynomial will be approximated from here")
print("Note: one of these lists is an iterable numpy array and the other is a simple python list")
n=int(input("Enter the degree of Lagrange polynomial to be calculated: "))
Px=sympify("0")
for i in range(n):
y=fx[i]
for j in range(n):
if(i!=j):
y=y*((x-datapoints[i])/(datapoints[j]-datapoints[i]))
Px=Px+y
print(Px)
xaxis=np.linspace(-2,2,100,endpoint=False)
f=[]
for val in xaxis:
f.append(calcfunc(val))
def calcpoly(a):
return Px.evalf(subs={x:a})
lP=[]
for val in xaxis:
lP.append(calcpoly(val))
f=np.array(f)
lP=np.array(lP)
plt.plot(xaxis,f,c='r')
#plt.plot(xaxis, lP)
plt.show()