数值分析第一次作业.docx
《数值分析第一次作业.docx》由会员分享，可在线阅读，更多相关《数值分析第一次作业.docx（18页珍藏版）》请在冰点文库上搜索。
数值分析第一次作业
数值分析第一次作业
%p55第1题主程序
%工院雷祺
%2015.10.29
clc;
clear;
symst;
x=[0.20.40.60.81];
y=[0.980.920.810.640.38];
%i=[0,1,10,11];
i=linspace（0,10,11）;
xx=0.2+0.08.*i;
%三次样条
[f3,yys]=mysplineN（x,y,0,0,xx）;
%牛顿
[P4,yyn]=myLanguage（x,y,xx）;
%画图
plot（x,y,'*',xx,yyn,'r-.',xx,yys,'--y','Linewidth',2）;
title（'牛顿插值、三次样条插值'）;
legend（'插值节点','牛顿插值','三次样条插值'）;
xlabel（'x'）;
ylabel（'y'）;
gridon;
boxon;
%显示插值多项式
fprintf（'牛顿插值多项式：
'）
P4=vpa（P4,6）
fprintf（'三次样条插值多项式：
'）
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
S=
0.446429*t^3-1.33929*t^2+0.45*t+0.925714
S=
-1.69643*t^3+2.51786*t^2-1.86429*t+1.38857
S=
2.58929*t^3-7.76786*t^2+6.36429*t-0.805714
%p55第2题主程序
%工院雷祺
%2015.10.29
clc;
clear;
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）;
plot（x,y,'*',xx,yy,xx,yyL,'r-.',xx,yys,'--y','Linewidth',2）;
title（'多项式插值、三次样条插值与原函数的比较'）;
legend（'插值节点','f（x）','多项式插值','三次样条插值'）;
xlabel（'x'）;
ylabel（'y'）;
gridon;
boxon;
subplot（1,2,2）;
plot（xx,abs（yy-yyL）,'r-.',xx,abs（yy-yys）,'--y','Linewidth',2）;
title（'多项式插值与三次样条插值的误差'）;
legend（'多项式插值','三次样条'）;
xlabel（'x'）;
ylabel（'E（y）'）;
gridon;
boxon;
运行结果：
输入n=10时
输入n=20时
%p55第3题主程序
%工院雷祺
%2015.10.29
clc;
clear;
symst;
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,'r-.',xx,yys,'--y','Linewidth',2）;
title（'拉格朗日插值、三次样条插值'）;
legend（'f（x）','L8（x）','S（x）'）;
xlabel（'x'）;
ylabel（'y'）;
gridon;
boxon;
subplot（2,2,3）;
plot（xx,abs（yy-yyL）,'r-.',xx,abs（yy-yys）,'--y','Linewidth',2）;
title（'拉格朗日插值与三次样条插值在[0,64]的误差'）;
legend（'拉格朗日','三次样条'）;
xlabel（'x'）;
ylabel（'E（y）'）;
gridon;
boxon;
subplot（2,2,4）;
plot（xx（1:
11）,abs（yy（1:
11）-yyL（1:
11））,'r-.',xx（1:
11）,abs（yy（1:
11）-yys（1:
11））,'--y','Linewidth',2）;
title（'拉格朗日插值与三次样条插值在[0,1]的误差'）;
legend（'拉格朗日','三次样条'）;
xlabel（'x'）;
ylabel（'E（y）'）;
gridon;
boxon;
fprintf（'牛顿插值多项式：
'）
L8=vpa（L8,6）
fprintf（'三次样条插值多项式：
'）
n=length（f3）;
fori=1:
n
S=poly2sym（f3（i,1:
4）,t）;
S=vpa（S,6）
end
运行结果：
牛顿插值多项式：
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
三次样条插值多项式：
S=
-0.0868256*t^3+1.08683*t
S=
0.0320461*t^3-0.356615*t^2+1.44344*t-0.118872
S=
-0.00273689*t^3+0.0607806*t^2-0.226141*t+2.10724
S=
0.000649371*t^3-0.0306484*t^2+0.596719*t-0.361343
S=
-0.000102098*t^3+0.00542215*t^2+0.0195906*t+2.71667
S=
0.00013415*t^3-0.0122964*t^2+0.462555*t-0.974697
S=
-0.000297962*t^3+0.0343716*t^2-1.2175*t+19.1859
S=
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
symst;
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;
end
f=0.0;
fori=1:
n
l=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;
f=f+l;%计算拉格朗日插值函数
simplify（f）;%化简
ifi==n
ifnargin==3
yc=subs（f,'t',x0）;%计算插值点的函数值
else
f=collect（f）;%将插值多项式展开
f=vpa（f,6）;%将插值多项式的系数化成6位精度的小数
end
end
end
function[f,yc]=Newton（x,y,x0）
%x，y为已知序列，x0待插值序列。
%f返回插值多项式，yc返回插值结果
%如果不输入待插值序列则输出插值多项式
%雷祺2015.10.20
symst;
if（length（x）==length（y））%输入序列检验
n=length（x）;
c（1:
n）=0.0;
else
disp（'x和y维数不等'）;
return;
end
f=y
（1）;
y1=0;
l=1;
fori=1:
n-1
forj=i+1:
n
y1（j）=（y（j）-y（i））/（x（j）-x（i））;%均差
end
c（i）=y1（i+1）;
l=l*（t-x（i））;
f=f+c（i）*l;%插值多项式
simplify（f）;%化简
y=y1;
if（i==n-1）
if（nargin>2）
yc=subs（f,'t',x0）;%计算插值序列的值
else
f=collect（f）;%插值多项式展开
f=vpa（f,6）;
end
end
end
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）
%雷祺2015.10.20
symst;
f=0.0;
f0=[];
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;
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;%lamda
（1）
A（n,n-1）=1;%u（n）
fori=2:
n-1
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
end
d
（1）=6/h
（1）*（ff
（1）-y_1）;
d（n）=6/h（n-1）*（y_n-ff（n-1））;
m=followup（A,d）;%用追赶法求解方程组
fori=1:
n-1
%计算插值多项式
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）;
f=collect（f）;%插值多项式展开
fh（i,1:
4）=sym2poly（f）;%提取系数
f=0.0;
end
xn=length（x0）;%需要插值的序列个数
forj=1:
xn%求插值结果
fori=1:
n-1
%若插值点超过范围，按最近的插值多项式插值
ifx0（j）（1）
f=poly2sym（fh（1,:
）,t）;
f0（j）=subs（f,'t',x0（j））;
f=0.0;
break;
end
ifx0（j）>x（n）
f=poly2sym（fh（n-1,:
）,t）;
f0（j）=subs（f,'t',x0（j））;
f=0.0;
break;
end
if（x（i）<=x0（j））&&（x（i+1）>=x0（j））
f=poly2sym（fh（i,:
）,t）;
f0（j）=subs（f,'t',x0（j））;
f=0.0;
break;
end
end
end
function[fh,f0]=mysplineN（x,y,y_1,y_n,x0）
%第二种边界条件下的三次样条插值，已知边界点二次导
%x，y为已知序列，y_1,y_n为边界的二阶导数，x0为待插值序列
%fh返回插值多项式，f0返回插值结果
%d
（1）=2*y_1;
%d（n）=2*y_n;
%雷祺2015.10.20
symst;
f=0.0;
f0=[];
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;
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）=0;%lamda
（1）
A（n,n-1）=0;%u（n）
fori=2:
n-1
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
end
d
（1）=2*y_1;
d（n）=2*y_n;
m=followup（A,d）;%用追赶法求解方程组
fori=1:
n-1
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）;
f=collect（f）;%插值多项式展开
fh（i,1:
4）=sym2poly（f）;%提取系数
f=0.0;
end
xn=length（x0）;%需要插值的序列个数
forj=1:
xn%求插值结果
fori=1:
n-1
%若插值点超过范围，按最近的插值多项式插值
ifx0（j）（1）
f=poly2sym（fh（1,:
）,t）;
f0（j）=subs（f,'t',x0（j））;
f=0.0;
break;
end
ifx0（j）>x（n）
f=poly2sym（fh（n-1,:
）,t）;
f0（j）=subs（f,'t',x0（j））;
f=0.0;
break;
end
if（x（i）<=x0（j））&&（x（i+1）>=x0（j））
f=poly2sym（fh（i,:
）,t）;
f0（j）=subs（f,'t',x0（j））;
f=0.0;
break;
end
end
end
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）;
end
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
fori=2:
（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）;
end
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）;
end
升级会员 