数值分析第一次作业.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