实习题答案p50,p239(数值分析(第五版))
一
1 求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序函数如下：
function [y,R,A,C,L]=newdscg(X,Y,x,M)
n=length(X); m=length(x);
for t=1:m
z=x(t); A=zeros(n,n);A(:,1)=Y';
s=0.0; p=1.0; q1=1.0; c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1)); end
q1=abs(q1*(z-X(j-1)));c1=c1*j;
end
C=A(n,n);q1=abs(q1*(z-X(n)));
for k=(n-1):-1:1
C=conv(C,poly(X(k)));
d=length(C);C(d)=C(d)+A(k,k);
end
y(k)= polyval(C, z);
end
R=M*q1/c1;L(k,:)=poly2sym(C);
首先将名为newdscg.m 的程序保存为M 文件，然后在MA TLAB 工作窗口输入程序 >> syms M,X=[0.2,0.4,0.6,0.8,1.0];
Y =[0.98,0.92,0.81,0.64,0.38];
x=0.9;
[y,R,A,C,P]=newdscg(X,Y,x,M)
运行后输出插值y=f(0.9)及其误差限公式R ，三阶牛顿插值多项式P 及其系数向量C ，差商的矩阵A 如下:
y =
0.5239
R =
(7*M)/800000
A =
0.9800 0 0 0 0
0.9200 -0.3000 0 0 0
0.8100 -0.5500 -0.6250 0 0
0.6400 -0.8500 -0.7500 -0.2083 0
0.3800 -1.3000 -1.1250 -0.6250 -0.5208
C =
-0.5208 0.8333 -1.1042 0.1917 0.9800
P =
- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/50
2:用拉格朗日插值法求其4次插值多项式的过程如下：
先编写一个M 文件来求其拉格朗日差值多项式，其程序如下：
function f=Language(x,y,x0)
syms t l ;
if(length(x)==length(y))
n=length(x);
else
disp('x 和y 的维数不相等！' );
return ; %检错
end
h=sym(0);for (i=1:n)
l=sym(y(i));
for (j=1:i-1)
l=l*(t-x(j))/(x(i)-x(j));
end ;
for (j=i+1:n)
l=l*(t-x(j))/(x(i)-x(j));
end ;
h=h+l;
end
simplify(h);
if (nargin==3)
f=subs(h,'t' ,x0); %计算插值点的函数值?
else
f=collect(h);
f=vpa(f,6);%将插值多项式的系数化成6位精度的小数?
end
然后将Language.m 保存为M 文件, 然后在MA TLAB 工作窗口输入程序:
>> x=[0.2 0.4 0.6 0.8 1.0];
>> y=[0.98 0.92 0.81 0.64 0.38];
>> f=Language(x,y,0.9)
其运行结果为：
f =
0.5239
然后在输入如下程序可得其4次差值多项式
>> f=Language(x,y)
其运行结果为：
f =
- 0.520833*t^4 + 0.833333*t^3 - 1.10417*t^2 + 0.191667*t + 0.98
最后输入
>> plot(x,y)
可得如下
图
二
1. 牛顿迭代法：
dfun.m:
function fun = dfun( x )%求原方程的导数f ′(x)
%UNTITLED4 Summary of this function goes here
% Detailed explanation goes here
syms t
f=diff(fun0(t));
fun=subs(f,t,x);
end
fun0.m
function fun = fun0(x)
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
fun=x^2-3*x+2-exp(x);%所要求的原方程
end
newton1.m
function [y ,k ] = newton1 (x0 ,n ,derta)%y 为迭
%代序列,k 为迭代次数,x0 初始位置,n 为最大迭代次数,derta 为
%容许误差
%UNTITLED5 Summary of this function goes here
% Detailed explanation goes here
k = 1 ;
y = [x0] ;
t = x0 - fun0(x0)./dfun(x0) ;
while abs (t - x0) >=derta
y=[y ,t ] ;
x0 = t ;
k = k +1 ;
t = x0 - fun0(x0)./dfun(x0) ;
if (k -1) > n
error ('n is full') ;
end
End
1. 其结果：>> [y ,k ] = newton1(1,40,0.000000001)
y =
1.0000 0.2689 0.2575 0.2575
k =4
2. 不动点迭代：
程序如下：
f=inline('(x^2+2-exp(x))/3');
x=feval(f,0.5);
disp(x);
Eps=1E-8;
i=1;
while 1
x0=x;
i=i+1;
x=feval(f,x);
disp(x);
if abs(x-x0)
disp(i)
break;
end
end
运行结果如下：
0.[**************]
0.[**************]
0.[**************]
0.[**************]
0.[**************]
0.[**************]
0.[**************]
0.[**************]
0.[**************]
9
3. 斯特芬森迭代法程序
format long;
f=inline('x-((x^2+2-exp(x))/3-x)^2/((((x^2+2-exp(x))/3)^2+2-exp((x^2+2-exp(x))/3))/3-2*(x^2+2-exp(x))/3+x)');
disp('x=');
x=feval(f,0.5);
disp(x);
Eps=1E-8;
i=1;
while 1
x0=x;
i=i+1;
x=feval(f,x);
disp(x);
if abs(x-x0)
disp(i)
break; end
end
运行结果
0.[**************]
0.[**************]
0.[**************]
0.[**************] 4