數(shù)值分析第一次作業(yè)

1、數(shù)值分析第一次作業(yè)%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('插值節(jié)點(diǎn)','牛頓插值','三次樣條插值');xlabel('x');ylabel('y');grid on;box on;%顯示插值多項(xiàng)式fprintf('牛頓插值多項(xiàng)式：')P4=vpa(P4,6)fprintf('三次樣條插值多項(xiàng)式：')n=length(f3);for i=1:n S=poly2sym(f3(i,1:4),t); S=vpa(S,6)end運(yùn)行結(jié)果：牛頓插值多項(xiàng)式： P4 = - 0.

3、520833*t4 + 0.833333*t3 - 1.10417*t2 + 0.191667*t + 0.98 三次樣條插值多項(xiàng)式： 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('請(qǐng)輸入n的大小=');% n=20;x=linspace(-1,1,n);y=1./(1+25.*x.2);xx=linspace(-1,1,200);%取200個(gè)點(diǎn)進(jìn)行插值yy=1./(1+25.*xx.2);% 三次樣條f,yys=mysplineN(x,y,0,0,xx);% yys=spline(x,y,xx);%系統(tǒng)三次樣條% 多項(xiàng)式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('多項(xiàng)式插值、三次樣條插值與原函數(shù)的比較');legend('插值節(jié)點(diǎn)','f(x)','多項(xiàng)式插值','三次樣條插值');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;多項(xiàng)式插值與三次樣條插值的誤差');legend('多項(xiàng)式插值','三次樣條');xlabel('x');ylabel('E(y)');grid on;box on;運(yùn)行結(jié)果：輸入n=10時(shí)輸入n=20時(shí)%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個(gè)插值點(diǎn)yy=sqrt(xx);% 三次樣條% yyy=spline(x,y,xx);y

7、_1=0;%0.5*0(-1/2);%初值給不了。y_n=0.5*64(-0.5);%邊界的一階導(dǎo)數(shù)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、項(xiàng)式：')L8=vpa(L8,6)fprintf('三次樣條插值多項(xiàng)式：')n=length(f3);for i=1:n S=poly2sym(f3(i,1:4),t); S=vpa(S,6)end 運(yùn)行結(jié)果：牛頓插值多項(xiàng)式： 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 三次樣條插值多項(xiàng)式： 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插值函數(shù)：function f,yc=myLanguage(x,y,x0)%x，y為已知序列，x0待插值序列。%f返回插值多項(xiàng)式，yc返回插值結(jié)果%如果不輸入待插值序列則輸出插值多項(xiàng)式%雷祺2015.10.20syms t;if(length(x)=length(y) n=length(x);else disp('x和y的維數(shù)不相等！'); 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);%計(jì)算拉格朗日基函數(shù) end; f=f+l;%計(jì)算拉格朗日插值函數(shù) simplify(f);%化簡(jiǎn) if i=n if nargin=3 yc=subs(f,'t',x0); %計(jì)算插值點(diǎn)的函數(shù)值 else f=collect(f);%將插值多項(xiàng)式展開 f=vpa(f,6);%將插值多項(xiàng)式的系數(shù)化成6位精度的小數(shù) end endendfunction f,yc=Newton(x,y,x0)%x，y為

14、已知序列，x0待插值序列。%f返回插值多項(xiàng)式，yc返回插值結(jié)果%如果不輸入待插值序列則輸出插值多項(xiàng)式%雷祺2015.10.20syms t;if(length(x)=length(y)%輸入序列檢驗(yàn) n=length(x); c(1:n)=0.0;else disp('x和y維數(shù)不等'); 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;%插值多項(xiàng)式 simplify(f);%化簡(jiǎn)

15、 y=y1; if(i=n-1) if(nargin>2) yc=subs(f,'t',x0);%計(jì)算插值序列的值 else f=collect(f);%插值多項(xiàng)式展開 f=vpa(f,6); end endendfunction fh,f0 = myspline(x,y,y_1,y_n,x0)%已知端點(diǎn)的一階導(dǎo)數(shù)的三次樣條插值%x，y為已知序列，y_1,y_n為邊界的一階導(dǎo)數(shù)，x0為待插值序列%fh返回插值多項(xiàng)式，f0返回插值結(jié)果% 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的維數(shù)不相等！'); return;end %維數(shù)檢查h=diff(x);%計(jì)算步長(zhǎng)ff=diff(y)./diff(x);%一階均差A(yù) = diag(2*ones(1,n);%求解m的系數(shù)矩陣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); %形成系數(shù)矩陣及向量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 %計(jì)算插值多項(xiàng)式 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);%插值多項(xiàng)式展開 fh(i,1:4)=sym2poly(f);%提取系數(shù) f=0.0;endxn=length(x0);%需要插值的序列個(gè)數(shù)for j=1:xn %求插值結(jié)果 for i=1:n-1 %若插值點(diǎn)超過(guò)范圍，按最近的插值多項(xiàng)式插值 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)%第二種邊界條件下的三次樣條插值，已知邊界點(diǎn)

20、二次導(dǎo)%x，y為已知序列，y_1,y_n為邊界的二階導(dǎo)數(shù)，x0為待插值序列%fh返回插值多項(xiàng)式，f0返回插值結(jié)果% 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的維數(shù)不相等！'); return;end %維數(shù)檢查h=diff(x);%計(jì)算步長(zhǎng)ff=diff(y)./diff(x);%一階均差A(yù) = diag(2*ones(1,n);%求解m的系數(shù)矩陣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); %這里我改過(guò)是不是寫反了 A(i, i+1) = lamda(i); %形成系數(shù)矩陣及向量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);%插值多項(xiàng)式展開 fh(i,1:4)=sym2poly(f);%提取系數(shù) f=0.0;endxn=length(x0);%需要插值的序列個(gè)數(shù)for j=1:xn %求插值結(jié)果 for i=1:n-1 %若插值