数值分析大作业三、四、五、六、七
大作业 三
1. 给定初值
0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用
程序.
解：Matlab 程序如下：
函数m 文件：fu.m
function Fu=fu(x)
Fu=x^3/3-x;
end
函数m 文件：dfu.m
function Fu=dfu(x)
Fu=x^2-1;
end
用Newton 法求根的通用程序Newton.m
clear;
x0=input('请输入初值x0:');
ep=input('请输入容许误差:');
flag=1;
while flag==1
x1=x0-fu(x0)/dfu(x0);
if abs(x1-x0) flag=0; end x0=x1; end fprintf('方程的一个近似解为：%f\n',x0); 寻找最大δ值的程序：Find.m clear eps=input('请输入搜索精度：'); ep=input('请输入容许误差：'); flag=1; k=0; x0=0; while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) flag1=0; end m=m+1; x0=x1; end if flag1==1||abs(x0)>=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +??=-= ?-???解：Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while n x=x0-f(x0)/subs(df(),x0); if abs(x-x0)>errorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)])（2）子函数非线性函数f function y=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end （3）子函数非线性函数的一阶导数df function y=df() syms x1 y=log((513+0.6651*x1)/(513-0.6651*x1))-x1/(1400*0.0918); y=diff(y); end 运行结果如下： 迭代次数: n=5 所求非零根: 正根x1=767.3861 负根x2=-767.3861 大作业 四 试编写MATLAB 函数实现Newton 插值，要求能输 出插值多项式. 对函数21()14f x x =+在区间[-5，5]上实现10次多项式插值. 分析：（1）输出插值多项式。 （2）在区间[-5,5]均匀插入99个节点，计算这些节点上函数f （x ）的近似值，并在同一图上画出原函数和插值多项式的图形。 （3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。 解：Matlab 程序代码如下： %此函数实现y=1/(1+4*x^2)的n 次Newton 插值，n 由调用函数时指定 %函数输出为插值结果的系数向量（行向量）和插值多项式 function [t y]=func5(n) x0=linspace(-5,5,n+1)'; y0=1./(1.+4.*x0.^2); b=zeros(1,n+1); for i=1:n+1 s=0; for j=1:i t=1; for k=1:i if k~=j t=(x0(j)-x0(k))*t; end ; end ; s=s+y0(j)/t; end; b(i)=s; end; t=linspace(0,0,n+1); for i=1:n s=linspace(0,0,n+1); s(n+1-i:n+1)=b(i+1).*poly(x0(1:i)); t=t+s; end; t(n+1)=t(n+1)+b(1); y=poly2sym(t); 10次插值运行结果： [b Y]=func5(10) b = Columns 1 through 4 -0.0000 0.0000 0.0027 -0.0000 Columns 5 through 8 -0.0514 -0.0000 0.3920 -0.0000 Columns 9 through 11 -1.1433 0.0000 1.0000 Y = - (10035*x^10)/2928 + x^9/616 + (256*x^8)/93425 - x^7/76 - (693*x^6)/1312 - (3*x^5)/ + (36624*x^4)/93425 - (5*x^3)/ - (32311*x^2)/70496 + (7*x)/ + 1