数值分析第一次作业文档格式.docx
《数值分析第一次作业文档格式.docx》由会员分享，可在线阅读，更多相关《数值分析第一次作业文档格式.docx（18页珍藏版）》请在冰点文库上搜索。
插值节点'
牛顿插值'
三次样条插值'
xlabel（'
x'
ylabel（'
y'
gridon;
boxon;
%显示插值多项式
fprintf（'
牛顿插值多项式：
'
）
P4=vpa（P4,6）
三次样条插值多项式：
n=length（f3）;
fori=1:
n
S=poly2sym（f3（i,1:
4）,t）;
S=vpa（S,6）
end
运行结果：
P4=
-0.520833*t^4+0.833333*t^3-1.10417*t^2+0.191667*t+0.98
S=
-1.33929*t^3+0.803571*t^2-0.407143*t+1.04
0.446429*t^3-1.33929*t^2+0.45*t+0.925714
-1.69643*t^3+2.51786*t^2-1.86429*t+1.38857
2.58929*t^3-7.76786*t^2+6.36429*t-0.805714
%p55第2题主程序
n=input（'
请输入n的大小='
%n=20;
x=linspace（-1,1,n）;
y=1./（1+25.*x.^2）;
xx=linspace（-1,1,200）;
%取200个点进行插值
yy=1./（1+25.*xx.^2）;
[f,yys]=mysplineN（x,y,0,0,xx）;
%yys=spline（x,y,xx）;
%系统三次样条
%多项式
[f,yyL]=myLanguage（x,y,xx）;
subplot（1,2,1）;
xx,yy,xx,yyL,'
多项式插值、三次样条插值与原函数的比较'
f（x）'
多项式插值'
subplot（1,2,2）;
plot（xx,abs（yy-yyL）,'
xx,abs（yy-yys）,'
多项式插值与三次样条插值的误差'
三次样条'
E（y）'
输入n=10时
输入n=20时
%p55第3题主程序
x=[01491625364964];
y=[012345678];
xx=linspace（0,64,641）;
%641个插值点
yy=sqrt（xx）;
%yyy=spline（x,y,xx）;
y_1=0;
%0.5*0^（-1/2）;
%初值给不了。
。
y_n=0.5*64^（-0.5）;
%边界的一阶导数
[f3,yys]=mysplineN（x,y,y_1,y_n,xx）;
%拉格朗日
[L8,yyL]=myLanguage（x,y,xx）;
subplot（'
position'
[0.13,0.55,0.8,0.3]）;
plot（xx,yy,xx,yyL,'
拉格朗日插值、三次样条插值'
L8（x）'
S（x）'
subplot（2,2,3）;
拉格朗日插值与三次样条插值在[0,64]的误差'
拉格朗日'
subplot（2,2,4）;
plot（xx（1:
11）,abs（yy（1:
11）-yyL（1:
11））,'
xx（1:
11）-yys（1:
拉格朗日插值与三次样条插值在[0,1]的误差'
L8=vpa（L8,6）
L8=
-3.28063e-10*t^8+6.71268e-8*t^7-0.00000542921*t^6+0.000222972*t^5-0.00498071*t^4+0.0604294*t^3-0.38141*t^2+1.32574*t
-0.0868256*t^3+1.08683*t
0.0320461*t^3-0.356615*t^2+1.44344*t-0.118872
-0.00273689*t^3+0.0607806*t^2-0.226141*t+2.10724
0.000649371*t^3-0.0306484*t^2+0.596719*t-0.361343
-0.000102098*t^3+0.00542215*t^2+0.0195906*t+2.71667
0.00013415*t^3-0.0122964*t^2+0.462555*t-0.974697
-0.000297962*t^3+0.0343716*t^2-1.2175*t+19.1859
0.000903973*t^3-0.142313*t^2+7.44004*t-122.221
插值函数：
function[f,yc]=myLanguage（x,y,x0）
%x，y为已知序列，x0待插值序列。
%f返回插值多项式，yc返回插值结果
%如果不输入待插值序列则输出插值多项式
%雷祺2015.10.20
if（length（x）==length（y））
n=length（x）;
else
disp（'
x和y的维数不相等！
return;
f=0.0;
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
end;
for（j=i+1:
n）
%计算拉格朗日基函数
f=f+l;
%计算拉格朗日插值函数
simplify（f）;
%化简
ifi==n
ifnargin==3
yc=subs（f,'
t'
x0）;
%计算插值点的函数值
else
f=collect（f）;
%将插值多项式展开
f=vpa（f,6）;
%将插值多项式的系数化成6位精度的小数
end
function[f,yc]=Newton（x,y,x0）
if（length（x）==length（y））%输入序列检验
c（1:
n）=0.0;
x和y维数不等'
f=y
（1）;
y1=0;
l=1;
n-1
forj=i+1:
y1（j）=（y（j）-y（i））/（x（j）-x（i））;
%均差
c（i）=y1（i+1）;
l=l*（t-x（i））;
f=f+c（i）*l;
%插值多项式
%化简
y=y1;
if（i==n-1）
if（nargin>
2）
%计算插值序列的值
%插值多项式展开
function[fh,f0]=myspline（x,y,y_1,y_n,x0）
%已知端点的一阶导数的三次样条插值
%x，y为已知序列，y_1,y_n为边界的一阶导数，x0为待插值序列
%fh返回插值多项式，f0返回插值结果
%d
（1）=6/h
（1）*（ff
（1）-y_1）;
%d（n）=6/h（n-1）*（y_n-ff（n-1））;
%A（1,2）=1;
%lamda
（1）
%A（n,n-1）=1;
%u（n）
f0=[];
if（length（x）==length（y））
n=length（x）;
end%维数检查
h=diff（x）;
%计算步长
ff=diff（y）./diff（x）;
%一阶均差
A=diag（2*ones（1,n））;
%求解m的系数矩阵
u=zeros（n-1,1）;
lamda=zeros（n-1,1）;
d=zeros（n,1）;
A（1,2）=1;
A（n,n-1）=1;
fori=2:
lamda（i）=h（i）/（h（i）+h（i-1））;
u（i-1）=h（i-1）/（h（i）+h（i-1））;
d（i）=6*（ff（i）-ff（i-1））/（h（i）+h（i-1））;
A（i,i-1）=u（i-1）;
A（i,i+1）=lamda（i）;
%形成系数矩阵及向量c
d
（1）=6/h
（1）*（ff
（1）-y_1）;
d（n）=6/h（n-1）*（y_n-ff（n-1））;
m=followup（A,d）;
%用追赶法求解方程组
%计算插值多项式
f=m（i）*（x（i+1）-t）^3/（6*h（i））+m（i+1）*（t-x（i））^3/（6*h（i））...
+（y（i）-m（i）*h（i）^2/6）*（x（i+1）-t）/h（i）+（y（i+1）-m（i+1）*h（i）^2/6）*（t-x（i））/h（i）;
fh（i,1:
4）=sym2poly（f）;
%提取系数
f=0.0;
xn=length（x0）;
%需要插值的序列个数
forj=1:
xn%求插值结果
fori=1:
%若插值点超过范围，按最近的插值多项式插值
ifx0（j）<
x
（1）
f=poly2sym（fh（1,:
）,t）;
f0（j）=subs（f,'
x0（j））;
break;
ifx0（j）>
x（n）
f=poly2sym（fh（n-1,:
if（x（i）<
=x0（j））&
&
（x（i+1）>
=x0（j））
f=poly2sym（fh（i,:
f0（j）=subs（f,'
function[fh,f0]=mysplineN（x,y,y_1,y_n,x0）
%第二种边界条件下的三次样条插值，已知边界点二次导
%x，y为已知序列，y_1,y_n为边界的二阶导数，x0为待插值序列
%d
（1）=2*y_1;
%d（n）=2*y_n;
%大小待定
A（1,2）=0;
A（n,n-1）=0;
%这里我改过是不是写反了
d
（1）=2*y_1;
d（n）=2*y_n;
f=0.0;
functionx=followup（A,f）%追赶法。
三次样条中用的
%AM=f
n=rank（A）;
for（i=1:
n-1）
a（i+1）=A（i+1,i）;
b（i）=A（i,i）;
c（i）=A（i,i+1）;
b（n）=A（n,n）;
%n=length（d）;
a
（1）=0;
%“追”的过程
L
（1）=b
（1）;
y
（1）=f
（1）/L
（1）;
u
（1）=c
（1）/L
（1）;
%beita1
（n-1）
L（i）=b（i）-a（i）*u（i-1）;
y（i）=（f（i）-y（i-1）*a（i））/L（i）;
u（i）=c（i）/L（i）;
L（n）=b（n）-a（n）*u（n-1）;
y（n）=（f（n）-y（n-1）*a（n））/L（n）;
%“赶”的过程
x（n）=y（n）;
fori=（n-1）:
-1:
1
x（i）=y（i）-u（i）*x（i+1）;
升级会员 