第二章 数值分析 上机作业 计科0801 邹庆华 学号 1081080126 第二章 第一题 牛顿插值函数 function y0=nutonpx(x0,y,xc) sizes=size(x0); n=sizes(2); f=zeros(n-1,n-1); for k=1:n-1 xd=x0(k+1)-x0(k); yd=y(k+1)-y(k); f(k,1)=yd/xd; end for k=2:n-1 for i=k:n-1 xd=x0(i+1)-x0(i-k+1); yd=f(i,k-1)-f(i-1,k-1); f(i,k)=yd/xd; end end syms x px; px=y(1); for k=2:n w=1; for i=1:k-1 w=w*(x-x0(i)); end px=px+w*f(k-1,k-1); end y0=subs(px,x,xc); 求差商函数 function f=junchabiao(x,y) %x=[0.40 0.55 0.65 0.80 0.90 1.05]; %y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25382]; sizes=size(x); n=sizes(2); f=zeros(n-1,n-1); for k=1:n-1 xd=x(k+1)-x(k); yd=y(k+1)-y(k); f(k,1)=yd/xd; end for k=2:n-1 for i=k:n-1 xd=x(i+1)-x(i-k+1); yd=f(i,k-1)-f(i-1,k-1); f(i,k)=yd/xd; end end 三次样条函数 function [y s]=spl(x0,y0,xt,v1,un,s0,sn) %x0=[0.25 0.30 0.39 0.45 0.53]; %y0=[0.50 0.5477 0.6245 0.6708 0.7280]; %v1=0; un=0; s0=0; sn=0; f0=s0; fn=sn; h=diff(x0); n=length(x0)-1; v=zeros(1,n); u=zeros(1,n); D=ones(1,n+1)*2; v(1)=v1; u(n)=un; for k=1:n-1 u(k)=h(k)/(h(k)+h(k+1)); v(k+1)=1-u(k); end d=ones(1,n+1); f=junchabiao(x0,y0); %d(1)=(f(1,1)-f0)*6/h(1); d(1)=f0; d(2:n)=f(2:n,2)*6; %d(n+1)=(fn-f(4,1))*6/h(n); d(n+1)=fn; A=diag(D); for k=1:n A(k,k+1)=v(k); A(k+1,k)=u(k); end b=d'; m=A\b; syms x s; for k=1:n s(k)=m(k)*(x0(k+1)-x).^3/(6*h(k)); s(k)=s(k)+ m(k+1)*(x-x0(k)).^3/(6*h(k)); s(k)=s(k)+(x0(k+1)-x)*(y0(k)-m(k).*h(k).^2/6)/h(k); s(k)=s(k)+(x-x0(k))*(y0(k+1)-m(k+1).*h(k).^2/6)/h(k); end %[A b m] for k=1:4 if xt>=x0(k) if xt<=x0(k+1) y=subs(s(k),x,xt); break; end end end if xt>x0(n+1) fprintf('ERROR ! x > %.2f\n',x0(n+1)); y=subs(s(n),x,xt); else if xt> ERROR ! x > 1.00 s = 1441/1400 - (75*(conj(x) - 1/5)^3)/56 - (69*conj(x))/280 (75*(conj(x) - 3/5)^3)/56 - (25*(conj(x) - 2/5)^3)/28 - (159*conj(x))/280 + 1621/1400 (25*(conj(x) - 4/5)^3)/28 - (145*(conj(x) - 3/5)^3)/56 - (219*conj(x))/280 + 1801/1400 Spline newton插值结果 如下图 第三章 第二题 拟合函数 function f_2(x0,y0,n) n=n+1; m=length(x0); hold on; grid; plot(x0,y0,'*'); syms x ; for k=1:n p(k)=x.^(k-1); end A=zeros(n); b=zeros(1,n); for k=1:n for i=1:n for j=1:m A(k,i)=A(k,i)+subs(p(k),x,x0(j))*subs(p(i),x,x0(j)); end end end for k=1:n for i=1:m b(k)=b(k)+subs(p(k),x,x0(i))*y0(i); end end s0=A\b'; s=p*s0; y=subs(s,x,x0); plot(x0,y,'r'); hold off; 输入 clear; x0=[0.0 0.1 0.2 0.3 0.5 0.8 1.0]; y0=[1.0 0.41 0.50 0.61 0.91 2.02 2.46]; n=3; f_2(x0,y0,n) 得下图： 输入 n=3, f_2(x0,y0,n) n=4; f_2(x0,y0,n) n=5 n=3,n=4 第五章 第二题 高斯列主消元函数 function y=f5_2(A,b) s=size(A); m=s(1) ; n=s(2) ; f=[A b]; if det(A)~=0 for i=1:m-1 k=i; max=f(i,i); for j=i+1:n if abs(f(j,i))>abs(max) max=f(j,i); k=j; end end if k~=i&f(j,i)~=0 p=f(i,:); f(i,:)=f(k,:); f(k,:)=p; end for j=i+1:n f(j,:)=f(j,:)-f(i,:)*f(j,i)/f(i,i); end end y=zeros(1,3); for k=1:n i=n+1-k; if i==n y(n)=f(n,m+1)/f(n,m); else a=f(i,n+1); for j=i+1:n a=a-f(i,j)*y(j); end y(i)=a/f(i,i); end end for i=1:n fprintf('\n x%d=% -.4f',i,y(i)); end else fprintf('\n ERROR! det(A)=0\n'); end (1) x1=-1.8535 x2= 0.8471 x3= 0.7392 A = 3.0100 6.0300 1.9900 1.2700 4.1600 -0.2300 0.9870 -4.8100 9.3400 b = 1 1 1 detA = 20.3992 x = -1.8535 0.8471 0.7392 condA = 45.4245 (2) clear; A=[ 3.00 6.03 1.99 1.27 4.16 -1.23 0.99 -4.81 9.34]; b=[1 1 1]'; x=f5_2(A,b); detA=det(A); condA=cond(A); A b detA x condA x1= 119.5273 x2=-47.1426 x3=-36.8403 A = 3.0000 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9900 -4.8100 9.3400 b = 1 1 1 detA = -0.4070 x = 119.5273 -47.1426 -36.8403 condA = 2.3028e+003 比较condA1 和condA2 可知 x1比x2 更精确 x2的波动比较大，属于病态。 x1的结果更理想 第七章 （1） 迭代法 1 clear; %x=x.^2-3*x+2-exp(x); syms x y; f=(x.^2+2-exp(x))/3; xk=1; e=1; n=0; while abs(e)>=10^-8 y=subs(f,x,xk); e=y-xk; xk=y; n=n+1; end xk n 结果 xk = 0.2575 n = 15 牛顿 迭代法 clear; %3x=x.^2-3*x+2-exp(x); syms x y; f=x.^2-3*x+2-exp(x); p=x-f/diff(f); xk=1; e=1; n=0; while abs(e)>=10^-8 y=subs(p,x,xk); e=y-xk; xk=y; n=n+1; end xk n 结果 xk = 0.2575 n = 4 15>4 显然 牛顿法更好 （2） clear; %x=(20-x.^3-2*x.^2)/10; syms x y; h=0.1; f=x.^3+2*x.^2+10*x-20; g=(x+h).^3+2*(x+h).^2+10*(x+h)-20; p=x-h*f/(g-f); xk=0; e=1; n=0; while abs(e)>=10^-8 y=subs(p,x,xk); e=y-xk; xk=y; n=n+1; end xk n 结果 xk = 1.3688 n = 8 牛顿 迭代法 clear; %x=(20-x.^3-2*x.^2)/10; syms x y; h=0.1; f=x.^3+2*x.^2+10*x-20; p=x-f/diff(f); xk=0; e=1; n=0; while abs(e)>=10^-8 y=subs(p,x,xk); e=y-xk; xk=y; n=n+1; end xk n 结果 xk = 1.3688 n = 6 因为 8>6 所以 牛顿迭代法的迭代速度更快 第九章 1 (1) 改进欧拉方法 function f=f9_1(x,y) f=1/x^2-y/x; function f91(x0,y0,h,n ) x=x0; y=y0; for k=1:n yp=y(k)+h*f9_1(x(k),y(k)); yc=y(k)+h*f9_1(x(k)+h,yp); y(k+1)=(yp+yc)/2; x(k+1)=x(k)+h; end fprintf('x0=%.4f y0=%.4f h=%.4f\n',x0,y0,h); fprintf('变量 x值 (改进欧拉方法)数值解 y值 \n'); for i=1:n/10:n+1 fprintf('%.4f %.4f\n',x(i), y(i)); end end >> f91(1,1,0.05 ,20) x0=1.0000 y0=1.0000 h=0.0500 变量 x值 (改进欧拉方法)数值解 y值 1.0000 1.0000 1.1000 0.9958 1.2000 0.9853 1.3000 0.9711 1.4000 0.9547 1.5000 0.9371 1.6000 0.9188 1.7000 0.9004 1.8000 0.8822 1.9000 0.8642 2.0000 0.8467 经典的四阶R-K方法 function f=f9_1(x,y) f=1/x^2-y/x; function f920(x0,y0,h,n ) x=x0; y=y0; for k=1:n k1=f9_1(x(k),y(k)); k2=f9_1(x(k)+h/2,y(k)+h/2*k1); k3=f9_1(x(k)+h/2,y(k)+h/2*k2); k4=f9_1(x(k)+h,y(k)+h*k3); y(k+1)=y(k)+h/6*(k1+k4+k2*2+k3*2); x(k+1)=x(k)+h; end fprintf('x0=%.4f y0=%.4f h=%.4f\n',x0,y0,h); fprintf('变量 x值 (R-K)数值解 y值 \n'); for i=1:n/10:n+1 fprintf('%.4f %.4f\n',x(i), y(i)); end end >> f920(1,1,0.1 ,10) x0=1.0000 y0=1.0000 h=0.1000 变量 x值 (R-K)数值解 y值 1.0000 1.0000 1.1000 0.9957 1.2000 0.9853 1.3000 0.9710 1.4000 0.9546 1.5000 0.9370 1.6000 0.9188 1.7000 0.9004 1.8000 0.8821 1.9000 0.8641 2.0000 0.8466 (2) 经典的四阶R-K方法 function f=f9_2(x,y) f=-50*y+50*x^2+2*x; function f92(x0,y0,h,n ) a=0:1/n:1; b=1/3*exp(-50*a)+a.^2; x=x0; y=y0; for k=1:n k1=f9_2(x(k),y(k)); k2=f9_2(x(k)+h/2,y(k)+h/2*k1); k3=f9_2(x(k)+h/2,y(k)+h/2*k2); k4=f9_2(x(k)+h,y(k)+h*k3); y(k+1)=y(k)+h/6*(k1+k4+k2*2+k3*2); x(k+1)=x(k)+h; end fprintf('x0=%.4f y0=%.4f h=%.4f\n',x0,y0,h); fprintf('变量 x值 原函数 y值 (R-K)数值解 y值 \n'); for i=1:n/10:n fprintf('%.4f %.4f %.4f\n',a(i), b(i), y(i)); end end >> f92(0,1/3,0.1 ,10) x0=0.0000 y0=0.3333 h=0.1000 变量 x值 原函数 y值 (R-K)数值解 y值 0.0000 0.3333 0.3333 0.1000 0.0122 4.6055 0.2000 0.0400 63.0625 0.3000 0.0900 864.0494 0.4000 0.1600 11843.6300 0.5000 0.2500 162354.5110 0.6000 0.3600 2225606.7134 0.7000 0.4900 30509354.2778 0.8000 0.6400 418232392.1740 0.9000 0.8100 5733269034.7809 >> f92(0,1/3,0.025 ,40) x0=0.0000 y0=0.3333 h=0.0250 变量 x值 原函数 y值 (R-K)数值解 y值 0.0000 0.3333 0.3333 0.1000 0.0122 0.0130 0.2000 0.0400 0.0401 0.3000 0.0900 0.0900 0.4000 0.1600 0.1600 0.5000 0.2500 0.2500 0.6000 0.3600 0.3600 0.7000 0.4900 0.4900 0.8000 0.6400 0.6400 0.9000 0.8100 0.8100 >> f92(0,1/3,0.01 ,100) x0=0.0000 y0=0.3333 h=0.0100 变量 x值 原函数 y值 (R-K)数值解 y值 0.0000 0.3333 0.3333 0.1000 0.0122 0.0123 0.2000 0.0400 0.0400 0.3000 0.0900 0.0900 0.4000 0.1600 0.1600 0.5000 0.2500 0.2500 0.6000 0.3600 0.3600 0.7000 0.4900 0.4900 0.8000 0.6400 0.6400 0.9000 0.8100 0.8100 原函数 clear; x=0:0.1:1; y=1/3*exp(-50*x)+x.^2; s=[x' y'] s = 0 0.3333 0.1000 0.0122 0.2000 0.0400 0.3000 0.0900 0.4000 0.1600 0.5000 0.2500 0.6000 0.3600 0.7000 0.4900 0.8000 0.6400 0.9000 0.8100 1.0000 1.0000 显然 当 h=0.1 的时候 结果 的误差很大 当 h=0.25 ,0.01时 比较准确 而 h=0.01 时，最接近原函数
本文档为【数值分析作业】，请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑，
图片更改请在作品中右键图片并更换，文字修改请直接点击文字进行修改，也可以新增和删除文档中的内容。
该文档来自用户分享，如有侵权行为请发邮件ishare@vip.sina.com联系网站客服，我们会及时删除。
[版权声明] 本站所有资料为用户分享产生，若发现您的权利被侵害，请联系客服邮件isharekefu@iask.cn，我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权，请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传，仅限个人学习分享使用，禁止用于任何广告和商用目的。
浙公网安备 33021202002483