درونیابی

کدنویسی الگوریتم‌های 

1. لاگرانژ
2. روش تفاضلات تقسیم شده‌ی نیوتن
3. درونیابی با نقاط متساوی الفاصله
1. روش تفاضلات پیشرو
2. روش تفاضلات پسرو
4. نکات

هر کدام از این برنامه‌ها، ورودی‌های مختص خود را دارند{بخش Input} و xi و fi معرف نقاط بکار برده شده برای درونیابی می‌باشند؛ و xm نقطه‌ایست که باید مقدار آن را پس از درونیابی بدست آوریم.

در این برنامه‌ها از ماژول Mpmath برای کنترل دقت محاسبات و ماژول‌های Numpy و Matplotlib برای رسم تابع درونیابی استفاده شده است.

 روش لاگرانژ

from mpmath import *
# mpmath Option
mp.dps = 5; mp.pretty = True
# Attention: mpmath have more priority from other libs.
from pylab import plot,show
from numpy import linspace

#------ Input ----------------
xi = [0, .12, .19, .32, .4, .51]
fi = [1, 1.3, 1.8, 2.2, 2.8, 3.2]
xm = .25
#------ Calc ----------------
def f(x):  
    global fi,xi
    return fi[xi.index(x)]

def L(j,x):
    global xi
    n = len(xi) 
    lst = [ mpf(x     - xi[i]) for i in range(n)]
    lmj = [ mpf(xi[j] - xi[i]) for i in range(n)]
    lst.pop(j)
    lmj.pop(j)
    return prd(lst)/prd(lmj)

def prd(l):
    p = 1
    for c in l:
        p *= c
    return p

n = len(xi)

px = lambda x: sum( fi[i] * L(i,x) \
                    for i in range(n) )

###------ Print & Plot -----------------
print('Result {:}'.format(px(xm)))

xpoints = []
ypoints = []
for x in linspace(xi[0],xi[-1],100):
    xpoints += [x]
    ypoints += [px(x)]
plot(xpoints, ypoints)
plot(xi,fi, "ro")
show()

• خروجی

Result 1.9908

• نمودار برنامه

---

 روش تفاضلات تقسیم شده‌ی نیوتن

from mpmath import *
# mpmath Option
mp.dps = 5; mp.pretty = True
# Attention: mpmath have more priority from other libs.
from pylab import plot,show
from numpy import linspace

#------ Input ----------------
xi = [0, .12, .19, .32, .4, .51]
fi = [1, 1.3, 1.8, 2.2, 2.8, 3.2]
xm = .25
#------ Calc ----------------
def f(x):  
    global fi,xi
    return fi[xi.index(x)]

def prd(l):
    p = 1
    for c in l:
        p *= c
    return p

def fr(l):
    global ls
    if len(l) == 1 :
        return f(l[0])
    rev = ( fr(l[1:]) - fr(l[:-1]) ) / ( l[-1] - l[0] )
    ls += [mpf(rev)]
    return rev

ls = []
fr(xi)
fc = [f(xi[0])] + ls[ -(len(xi) - 1) :]
n = len(fc)

[print(" {:<{width}}".format(fc[i],width=16)) \
 for i in range(n)]

px = lambda x: sum( fc[i] * prd( [x - xi[j] for j in range(i)] ) \
                    for i in range(n) )

###------ Print & Plot -----------------
print('\nResult {:}'.format(px(xm)))

xpoints = []
ypoints = []
for x in linspace(xi[0],xi[-1],100):
    xpoints += [x]
    ypoints += [px(x)]
plot(xpoints, ypoints)
plot(xi,fi, "ro")
show()

• خروجی

 1               
 2.5             
 24.436          
 -139.89         
 719.3           
 -2804.0         

Result 1.9908

• نمودار برنامه

---

---

حالت خاص: درونیابی با نقاط متساوی الفاصله

 روش تفاضلات پیشرو

from mpmath import *
# mpmath Option
mp.dps = 5; mp.pretty = True
# Attention: mpmath have more priority from other libs.
from pylab import plot,show
from numpy import linspace

#------ Input ----------------
xi = [0, .1, .2, .3, .4, .5]
fi = [1, 1.3, 1.8, 2.2, 2.8, 3.2]
xm = .25
#------ Calc ----------------
def f(x):  
    global fi,xi
    return fi[xi.index(x)]

def c(s, k):
    s = mpf(s)
    if s-k < 0 and s-int(s) == 0: return 0
    return fac(s) / ( fac(k) * fac(s - k) )

def fr(l):
    global ls
    if len(l) == 1 :
        return f(l[0])
    rev = fr(l[1:]) - fr(l[:-1])
    ls += [mpf(rev)]
    return rev

ls = []
fr(xi)
fc = [f(xi[0])] + ls[ -(len(xi) - 1) :]
n = len(fc)
h = abs(xi[1] - xi[0])

[print(" {:<{width}}".format(fc[i],width=16)) \
 for i in range(n)]

px = lambda x: sum( fc[i] * c((x - xi[0]) / h, i) \
                    for i in range(n) )

###------ Print & Plot -----------------
print('\nResult {:}'.format(px(xm)))

xpoints = []
ypoints = []
for x in linspace(xi[0],xi[-1],100):
    xpoints += [x]
    ypoints += [px(x)]
plot(xpoints, ypoints)
plot(xi,fi, "ro")
show()

• خروجی

 1               
 0.3             
 0.2             
 -0.3            
 0.6             
 -1.3            

Result 1.9926

• نمودار برنامه

---

 روش تفاضلات پسرو

from mpmath import *
# mpmath Option
mp.dps = 5; mp.pretty = True
# Attention: mpmath have more priority from other libs.
from pylab import plot,show
from numpy import linspace

#------ Input ----------------
xi = [1.25, 1.26, 1.27, 1.28, 1.29]
fi = [1.3914, 1.3770, 1.3478, 1.3046, 1.2477]
xm = 1.265
#------ Calc ----------------
def f(x):  
    global fi,xi
    return fi[xi.index(x)]

def c(s, k):
    s = mpf(s)
    if s-k < 0 and s-int(s) == 0: return 0
    return fac(s) / ( fac(k) * fac(s - k) )

def fr(l):
    global ls
    if len(l) == 1 :
        return f(l[0])
    rev = -(fr(l[:-1]) - fr(l[1:]))
    ls += [mpf(rev)]
    return rev

ls = []
fr(xi)
fc = [f(xi[-1])] + ls[ -(len(xi) - 1) :]
n = len(fc)
h = abs(xi[1] - xi[0])

[print(" {:<{width}}".format(fc[i],width=16)) \
 for i in range(n)]

px = lambda x: sum((-1)**(i) * fc[i] * c(-(x - xi[-1]) / h, i) \
                   for i in range(n) )

###------ Print & Plot -----------------
print('\nResult {:}'.format(px(xm)))

xpoints = []
ypoints = []
for x in linspace(xi[0],xi[-1],100):
    xpoints += [x]
    ypoints += [px(x)]
plot(xpoints, ypoints)
plot(xi,fi, "ro")
show()

• خروجی

 1.2477          
 -0.0569         
 -0.0137         
 0.0003          
 -0.0005         

Result 1.3642

• نمودار برنامه

---

 نکات:

• همانطور که مشاهده می‌شود؛ الگوریتم لاگرانژ از سادگی نسبتاً خوبی برخوردار است اما این سادگی با پیچیدگی زمانی بالایی همراه است. ضمناً داده‌های ورودی روش لاگرانژ و تفاضلات تقسیم شده یکی است؛ خروجی آنها نیز یکی است.
• الگوریتم‌های پیشرو و پسرو در خطوط ۲۶ ، ۳۲ و ۳۹ تفاوت دارند!
• الگوریتم‌های تقسیم شده و پیشرو در خطوط نقره‌ای رنگ با یکدیگر تفاوت دارند!
• در هر سه الگوریتم تابع fr به صورت بازگشتی عمل می‌کند
• ارزش ریاضی خط ۲۶ در کدها تفاضل پیشرو و پسرو یکسان است. و این جابجایی صرفاً به دلیل ایجاد ترتیب معکوس در فراخوانی بازگشتی انجام شده است تا در نهایت در خط ۳۲ دستور «

ls[ -(len(xi) - 1) :]

» داده‌های مطلوب روش را فراهم کند.