数值分析大作业三四五六七完整版

数值分析大作业三四五

六七

HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】

大

作

业 三

1.

给定初值

0x 及容许误差

，编制牛顿法解方程f (x )=0的通用

程序.

解：Matlab 程序如下：

函数m 文件：

function Fu=fu(x) Fu=x^3/3-x; end

函数m 文件：

function Fu=dfu(x) Fu=x^2-1; end

用Newton 法求根的通用程序 clear;

x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1;

while flag==1

x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

fprintf('方程的一个近似解为：%f\n',x0); 寻找最大δ值的程序： 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)

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+*x)/*x))-x/(1400*; end

（3）子函数 非线性函数的一阶导数df function y=df() syms x1

y=log((513+*x1)/*x1))-x1/(1400*; y=diff(y);

运行结果如下：

迭代次数: n=5

所求非零根: 正根x1= 负根x2=

大作业 四

试编写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)

Columns 1 through 4

Columns 5 through 8

Columns 9 through 11

Y =

b为插值多项式系数向量，Y为插值多项式。插值近似值：

x1=linspace(-5,5,101);

x=x1(2:100);

y=polyval(b,x)

y =

Columns 1 through 12

Columns 13 through 24

Columns 25 through 36

Columns 37 through 48

Columns 49 through 60

Columns 61 through 72

Columns 73 through 84

Columns 85 through 96

Columns 97 through 99

绘制原函数和拟合多项式的图形代码：

plot(x,1./(1+4.*x.^2))

hold all

plot(x,y,'r')

xlabel('X')

ylabel('Y')

title('Runge现象')

gtext('原函数')

gtext('十次牛顿插值多项式')

绘制结果：

误差计数并绘制误差图：

hold off

ey=1./(1+4.*x.^2)-y

ey =