1. function newton_divided_difference( f, min_x, max_x, n ) 2. %UNTITLED 使用牛顿差商法插值 3. % f为原函数,min_x max_x定义了区间,n为多项式次数 4. 5. %步骤1:求x,f,f[1],f[2]...f[n] 6. x = min_x : (max_x-min_x)/n : max_x; 7. for i = 1:1:n+1 8. y(i) = f(x(i)); 9. end 10. a=[x',y']; 11. for i = 1:1:n 12. for j = 1:1:(n-i+1) 13. a(j,i+2) = (a(j+1,i+1)-a(j,i+1))/(a(j+i,1)-a(j,1)); 14. end 15. end 16. 17. %步骤2:画原始函数图像 18. figure(1); 19. STEP = 0.001; 20. x = min_x : STEP : max_x; 21. for i = 1:1:((max_x-min_x)/STEP+1) 22. y_1(i) = f(x(i)); 23. end 24. plot(x,y_1,'r') 25. 26. %步骤3:画插值函数图像 27. hold on 28. for i = 1:1:((max_x-min_x)/STEP+1) 29. y_2(i) = a(1,n+2); 30. for j = n:-1:1 31. y_2(i) = y_2(i) * (x(i) - a(j,1)) + a(1,j+1); 32. end 33. end 34. plot(x,y_2,'b') 35. hold off 36. 37. %步骤4:画误差图像 38. figure(2); 39. y_error = y_1 - y_2; 40. plot(x,y_error) 41. end