数值分析第一次作业

1、数值分析第一次作业%p55第1题主程序%工院雷祺%2015.10.29clc;clear;syms t;x=0.2 0.4 0.6 0.8 1;y=0.98 0.92 0.81 0.64 0.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)

2、;title('牛顿插值、三次样条插值');legend('插值节点','牛顿插值','三次样条插值');xlabel('x');ylabel('y');grid on;box on;%显示插值多项式fprintf('牛顿插值多项式：')P4=vpa(P4,6)fprintf('三次样条插值多项式：')n=length(f3);for i=1:n S=poly2sym(f3(i,1:4),t); S=vpa(S,6)end运行结果：牛顿插值多项式： P4 = - 0.

3、520833*t4 + 0.833333*t3 - 1.10417*t2 + 0.191667*t + 0.98 三次样条插值多项式： S = - 1.33929*t3 + 0.803571*t2 - 0.407143*t + 1.04 S = 0.446429*t3 - 1.33929*t2 + 0.45*t + 0.925714 S = - 1.69643*t3 + 2.51786*t2 - 1.86429*t + 1.38857 S = 2.58929*t3 - 7.76786*t2 + 6.36429*t - 0.805714%p55第2题主程序%工院雷祺%2015.10.29clc;

4、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,'-

5、y','Linewidth',2);title('多项式插值、三次样条插值与原函数的比较');legend('插值节点','f(x)','多项式插值','三次样条插值');xlabel('x');ylabel('y');grid on;box on;subplot(1,2,2);plot(xx,abs(yy-yyL),'r-.',xx,abs(yy-yys),'-y','Linewidth',2);title(&

6、#39;多项式插值与三次样条插值的误差');legend('多项式插值','三次样条');xlabel('x');ylabel('E(y)');grid on;box on;运行结果：输入n=10时输入n=20时%p55第3题主程序%工院雷祺%2015.10.29clc;clear;syms t;x=0 1 4 9 16 25 36 49 64;y=0 1 2 3 4 5 6 7 8;xx=linspace(0,64,641);%641个插值点yy=sqrt(xx);% 三次样条% yyy=spline(x,y,xx);y

7、_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(

8、x)','S(x)');xlabel('x');ylabel('y');grid on;box on;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)');g

9、rid on;box on;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)');grid on;box on;fprintf('牛顿插值多

10、项式：')L8=vpa(L8,6)fprintf('三次样条插值多项式：')n=length(f3);for i=1:n S=poly2sym(f3(i,1:4),t); S=vpa(S,6)end 运行结果：牛顿插值多项式： L8 = - 3.28063e-10*t8 + 6.71268e-8*t7 - 0.00000542921*t6 + 0.000222972*t5 - 0.00498071*t4 + 0.0604294*t3 - 0.38141*t2 + 1.32574*t 三次样条插值多项式： S = - 0.0868256*t3 + 1.08683*t S

11、= 0.0320461*t3 - 0.356615*t2 + 1.44344*t - 0.118872 S = - 0.00273689*t3 + 0.0607806*t2 - 0.226141*t + 2.10724 S = 0.000649371*t3 - 0.0306484*t2 + 0.596719*t - 0.361343 S = - 0.000102098*t3 + 0.00542215*t2 + 0.0195906*t + 2.71667 S = 0.00013415*t3 - 0.0122964*t2 + 0.462555*t - 0.974697 S = - 0.000297

12、962*t3 + 0.0343716*t2 - 1.2175*t + 19.1859 S = 0.000903973*t3 - 0.142313*t2 + 7.44004*t - 122.221插值函数：function f,yc=myLanguage(x,y,x0)%x，y为已知序列，x0待插值序列。%f返回插值多项式，yc返回插值结果%如果不输入待插值序列则输出插值多项式%雷祺2015.10.20syms t;if(length(x)=length(y) n=length(x);else disp('x和y的维数不相等！'); return;endf = 0.0;for i

13、=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);%化简 if i=n if nargin=3 yc=subs(f,'t',x0); %计算插值点的函数值 else f=collect(f);%将插值多项式展开 f=vpa(f,6);%将插值多项式的系数化成6位精度的小数 end endendfunction f,yc=Newton(x,y,x0)%x，y为

14、已知序列，x0待插值序列。%f返回插值多项式，yc返回插值结果%如果不输入待插值序列则输出插值多项式%雷祺2015.10.20syms t;if(length(x)=length(y)%输入序列检验 n=length(x); c(1:n)=0.0;else disp('x和y维数不等'); return;endf=y(1);y1=0;l=1;for i=1:n-1 for j=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);%化简

15、 y=y1; if(i=n-1) if(nargin>2) yc=subs(f,'t',x0);%计算插值序列的值 else f=collect(f);%插值多项式展开 f=vpa(f,6); end endendfunction 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;%

16、lamda(1)% A(n,n-1)=1;%u(n)%雷祺2015.10.20syms t;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;

17、%u(n)for i=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); %形成系数矩阵及向量cendd(1)=6/h(1)*(ff(1)-y_1);d(n)=6/h(n-1)*(y_n-ff(n-1);m = followup(A,d);%用追赶法求解方程组for i=1:n-1 %计算插值多项式 f=m(i)*(x(i+1)-t)3/(6*h(i)+m(i+1)*(

18、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;endxn=length(x0);%需要插值的序列个数for j=1:xn %求插值结果 for i=1:n-1 %若插值点超过范围，按最近的插值多项式插值 if x0(j)<x(1) f=poly2sym(fh(1,:),t); f0(j)=subs(f,'t',x0(j); f

19、=0.0; break; end if x0(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 endendfunction fh,f0 = mysplineN(x,y,y_1,y_n,x0)%第二种边界条件下的三次样条插值，已知边界点

20、二次导%x，y为已知序列，y_1,y_n为边界的二阶导数，x0为待插值序列%fh返回插值多项式，f0返回插值结果% d(1) = 2*y_1;% d(n) = 2*y_n;%雷祺2015.10.20syms t;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);%大小待定

21、lamda = zeros(n-1,1);d = zeros(n,1);A(1,2)=0;%lamda(1)A(n,n-1)=0;%u(n)for i=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); %形成系数矩阵及向量cendd(1) = 2*y_1;d(n) = 2*y_n;m = followup(A,d);%用追赶法求解方程组for i=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;endxn=length(x0);%需要插值的序列个数for j=1:xn %求插值结果 for i=1:n-1 %若插值