p0=p; b=0 ; %工作变量
W=poly(X); %W为以X为根的多项式
dW=polyder(W); %dW 为对多项式W 求导后的多项式系数 z=polyval(dW,X); %z 为以 dW 为系数的多项式对X 的值 A=[1 1]; r =A ; %A,r 为长度为 2的(1,2)向量 for i=1:n %计算循环开始
A=[1,-X (i)]; %A 为一次多项式x-x(i)系数 [p0,r]=deconv(W,A); %进行多项式除法W/A,p0 为商 b=Y (i)/z(i);
p0=b.*p0 ; %p0为累加项 p=p+p0 ;
end %循环结束 disp(toc) %显示所用时间
(2) 例题
%建立Matlab命令文件 LagSin.m
%用拉格朗日插值法拟合 [0,2*pi]上的sin函数 程序命令行：
x0=linspace(0,2*pi,10);
y0=sin(x0); %拟合数据采样
p=Lag_polyfit(x0,y0); %调用Lag_polyfit 计算拟合多项式 x=0 :0.2:2*pi;
y1=sin(x); %sin 函数值 y2=polyval(p,x); %拟合多项式值
plot(x,y1,x,y2,'r') %sin 函数(蓝色)与多项式(红色)效果图 命令窗口中输入 lagsin观察图形输出窗口
>>lagsin %函数(文件)名不区分大小写
3．牛顿插值法
(1) 均差表计算及算法实现 function [p]=Newton_Polyfit(X,Y) %Matlab函数文件 Newton_polyfit.m %牛顿插值法多项式拟合
16
%[p]=Newton_Polyfit(X,Y) % X 拟合自变量 % Y 拟合函数值
% p 所得的拟合多项式系数 if size(X)~=size(Y)
error('变量不匹配?’);%如果要拟合的函数值与自变量维数不一样,则退出报错 end
format long g %设置最合适的数字格式 r=size(X); n=r(2); %n 为要拟合的数据长度 M=ones(n,n); %均差表工作矩阵
M(:,1)=Y'; %设置均差表第一列为函数值 for i=2:n
for j=i:n %高阶均差推算 M(j,i)=(M(j,i-1)-M(j-1,i-1))/ (X (j)-X (j-i+1)); end end
M %显示均差表 p0=[zeros(1,n-1) M(1,1)];p=p0 ; %拟合多项式存储向量及初值 for i=1:n-1
p1=M(i+1,i+1).*poly(X (1:i)); %对应阶均差的多项式 p0=[zeros(1,n-i-1) p1];
p=p+p0 ; %对应项累加 end (2) 例题
%Matlab 命令文件 NewtonSin.m %用牛顿插值法拟合 [0,2*pi]上的 sin函数 x0=linspace(0,2*pi,10);
y0=sin(x0); %拟合数据采样
p=Newton_polyfit(x0,y0); %调用Newton_polyfit 计算拟合多项式 x=0 :0.2:2*pi;
y1=sin(x); %sin 函数值 y2=polyval(p,x); %拟合多项式值
plot(x,y1,x,y2,'r') %sin 函数(蓝色)与多项式(红色)效果图 命令窗口中输入 newtonsin观察图形输出窗口
17
>>newtonsin
数值分析实验四
第五章 数值积分
1. 实验目的
1) 熟悉Matlab矩阵操作
2) 用Matlab实现数值积分Cores,Romberg及复化Gauss算法 3) 学会数值积分有关应用
2. 实验内容
(1) Cores积分
在使用牛顿-柯斯特公式时,通过提高阶为提高精度的途径并不能取得满意的效果.通常采用复化积分法.以复化柯斯特公式为例: 将每个区间 4 等分, 依此用柯斯特公式 function [val]=Cores_int(Y,u) % [val]=Cores_int(Y,u) %Y积分函数数值 % u区间间隔 % val返回积分值
r=size(Y); n=r(2); v=n; h=4*u; m=mod(n-1,4);
if m==0 %如果取值点恰巧满足复化区间数 t=(n-1)/4; x=ones(t,1);
val=(h/90)*(7*Y(1)+32*(Y(2:4:n-3)*x)+12*(Y(3:4:n-2)*x)+32*(Y(4:4:n-1)*x)+ ….
+14*(Y(5:4:n-4)*ones(size((5:4:n-4)')))+7*Y(n)); elseif m==1 n=n-1;
t=(n-1)/4; x=ones(t,1);
val=(h/90)*(7*Y(1)+32*(Y(2:4:n-3)*x)+12*(Y(3:4:n-2)*x)+32*(Y(4:4:n-1)*x)+ …. +14*(Y(5:4:n-4)*ones(size((5:4:n-4)')))+7*Y(n)); n=v;
val=val+(u/2)*(Y(n-1)+Y(n)); %多余的区间用梯形公式 elseif m==2 n=n-2;
18
t=(n-1)/4; x=ones(t,1);
val=(h/90)*(7*Y(1)+32*(Y(2:4:n-3)*x)+ 12*(Y(3:4:n-2)*x)+32*(Y(4:4:n- 1)*x)+14*(Y(5:4:n-4)*ones(size((5:4:n-4 )')))+ 7*Y(n)); n=v;
val=val+(u/3)*(Y(n-2)+4*Y(n-1)+Y(n)); %Simpson公式 elseif m==3 n=n-3;
t=(n-1)/4; x=ones(t,1);
val=(h/90)*(7*Y(1)+32*(Y(2:4:n-3)*x)+...
12*(Y(3:4:n-2)*x)+32*(Y(4:4:n-1)*x)+14*(Y(5:4:n-4)*ones(size((5:4:n-4)')))+7*Y(n)); n=v;
val=val+(3*u/8)*(Y(n-3)+3*Y(n-2)+3*Y(n-1)+Y(n)); %3阶Newton-Cores公式 end
例: 对于函数 f(x)=sinx/x,利用上述算法计算[0,1]区间上积分. >>x=0:1/8:1; >>x=x+eps; >>y=sin(x)./x;
>>[val]=Cores_int(y,1/8) val =
0.946083069350917 (2) Romberg算法:
1hn?1 复化梯形积分值递推公式: T2n?Tn??f22k?0km??x?k?1? ?2?(k?1,2,3)
4m1k?1k?TmRichardson外推加速法递推公式: T?mTm?1?1m4?14?1
function [val,M]=Romberg(f,ab,dalt) % Romberg算法计算积分 % [val,M]=Romberg(f,ab,dalt) %val返回积分值 %M T数表
% f被积函数,函数文件名 % ab积分区间 % dalt精度
19
if nargin<3
dalt=1e-7; %设置默认精度 10 end
i=1; h=ab(2)-ab(1); t=0; j=0; T(1,1)=(h/2)*(feval(f,ab)*ones(2,1));
while i<50 %最大跌代次数为 50次 i=i+1; h=h/2;
T=[T zeros(i-1,1);zeros(1,i)]; x=ab(1)+h:2*h:ab(2)-h;
y=feval(f,x)*ones(size(x')); % feval(f,x)求得 f在x点值 T(i,1)=T(i-1,1)/2+h*y; % 细化区间,求得梯形值 for t=2:i j=t-1;
T(i,t)=(4^j*T(i,j)-T(i-1,j))/(4^j-1); %外推加速 end
if abs(T(i,i)-T(i-1,i-1))<=dalt %控制精度 break end end
if nargout==2
M=T; %按需求返回 T数表 end
val=T(i,i); %积分值
例1 : 利用上述算法计算积分: 步骤:
1) 新建函数fun1.m function y=fun1(x) y=sqrt(1-sin(x).^2);
2) >> [val,M]=Romberg('fun1',[0 pi/2],1e-10)
20
??/201?sin2?d?