2019年1月17日 星期四

範例 1-9 依牛頓內插法的插除符號理則,建立差除表。

範例 1-9 已知6點座標如下:

n=5
0.0  -6.0
0.1  -5.89483
0.3  -5.65014
0.6  -5.17788
1.0  -4.28172
1.1  -3.99583

請依牛頓內插法的插除符號理則,建立差除表。

#C:Python3.7 python.exe (32bit)
# -*- coding:utf-8 -*-

'''
/* ex1-9.c is for developing the divided-defference
 * table for Newton Interpolation polynomial.
 */

n=5
0.0  -6.0
0.1  -5.89483
0.3  -5.65014
0.6  -5.17788
1.0  -4.28172
1.1  -3.99583

'''
from array import *
xa=1.5

x=[0.0 ,0.1 ,0.3 ,0.6 ,1.0 ,1.1]
f=[[-6.0 ,0.0,0.0,0.0,0.0,0.0,0.0] ,
   [-5.89483,0.0,0.0,0.0,0.0,0.0] ,
   [-5.65014,0.0,0.0,0.0,0.0] ,
   [-5.17788,0.0,0.0,0.0] ,
   [-4.28172,0.0,0.0] ,
   [-3.99583,0.0]]

# 留意 f Array的設定

result=0.0
n=5   # 0...5=6

print('x=',x)
print('f=',f)

print("\n\t\tDivided Difference Table:\n")
print(" ============================================\n")
 
for j in range (1,n+1):
    for i in range (0,n-j+1):
        f[i][j] =  (  (f[i+1][j-1]) - (f[i][j-1]) )  /  ( x[i+j] - x[i] )
     
print("i\tx(i)\tf(i)\tf(i,i+1)\tf(i,i+1,i+2)\tf(i,i+1,i+2,i+3) ......\n")
for i in range (0,n+1):
    print('%d ,\t%2.4f' %(i,x[i]),end='')
    for j in range(0,n-i+1):
        print('\t%8.5f' %f[i][j],end='')
    print("\n")




輸出畫面
x= [0.0, 0.1, 0.3, 0.6, 1.0, 1.1]
f= [[-6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [-5.89483, 0.0, 0.0, 0.0, 0.0, 0.0], [-5.65014, 0.0, 0.0, 0.0, 0.0], [-5.17788, 0.0, 0.0, 0.0], [-4.28172, 0.0, 0.0], [-3.99583, 0.0]]

Divided Difference Table:

 ============================================

i x(i) f(i) f(i,i+1) f(i,i+1,i+2) f(i,i+1,i+2,i+3) ......

0 , 0.0000 -6.00000 1.05170 0.57250 0.21500 0.06302 0.01416

1 , 0.1000 -5.89483 1.22345 0.70150 0.27802 0.07859

2 , 0.3000 -5.65014 1.57420 0.95171 0.35661

3 , 0.6000 -5.17788 2.24040 1.23700

4 , 1.0000 -4.28172 2.85890

5 , 1.1000 -3.99583

>>> 

=======================================================
原來C語言程式

/* ex1-9.c is for developing the divided-defference
 * table for Newton Interpolation polynomial.
 */
#include <stdio.h>
void main()
{
   int i,j,n;
   double f[40][40],x[40];
   scanf("n=%d",&n);
   for(i=0;i<=n;i++)
   {
      scanf("%lf %lf ",&x[i],&f[i][0]);
   }
   printf("\n Divided Difference Table:\n");
   printf(" =========================\n");
   for(j=1;j<=n;j++)
   {
      for(i=0;i<=n-j;i++)
      {
f[i][j]=(f[i+1][j-1]-f[i][j-1])/(x[i+j]-x[i]);
      }
   }
   printf("i   x(i)     f(i)    f(i,i+1) f(i,i+1.i+2),  ......\n");
   for(i=0;i<=n;i++)
   {
      printf("%d  %8.5lf ",i,x[i]);
      for(j=0;j<=n-i;j++)
      {
printf("%8.5lf ",f[i][j]);
      }
      printf("\n");
   }
   return;
}


沒有留言:

張貼留言

Messaging API作為替代方案

  LINE超好用功能要沒了!LINE Notify明年3月底終止服務,有什麼替代方案? LINE Notify將於2025年3月31日結束服務,官方建議改用Messaging API作為替代方案。 //CHANNEL_ACCESS_TOKEN = 'Messaging ...