数值分析编程实例-1.doc

数值计算方法上机题目1、利用Steffensen方法和Muller(抛物线)方法计算下列方程在区间上的实根，要求精度为，并比较迭代次数。Steffensen法:function A=steffensen(f,x,ep)%f为迭代函数%x为初始值%ep为精度A=;A(1,1)=x;for i=1:1000 y=feval(f,A(i,1); z=feval(f,y); A(i,2)=A(i,1)-(y-A(i,1).2/(z-2*y+A(i,1); if abs(A(i,2)-A(i,1)ep break end A(i+1,1)=A(i,2);end执行:f=inline(log(6-2*cos(x)-2(-x);y=1;ep=10.-5;steffensen(f,y,ep)结果:ans = 1.0000 1.8682 1.8682 1.8295 1.8295 1.8294 1.8294 1.8294Muller法:functionA=muller(f,x0,x1,x2,ep)%f为原函数%x0为输入初值最小的那个%x1为输入初值最大的那个%x2为输入初值中间的那个%ep为输入精度A=;A(1,1)=x0;A(1,2)=x1;A(1,3)=x2;for i=1:1000 y=(A(i,3)-A(i,2)/(A(i,2)-A(i,1); x=1+y; a=feval(f,A(i,1)*y.2-feval(f,A(i,2)*y*x+feval(f,A(i,3)*y; b=feval(f,A(i,1)*y.2-feval(f,A(i,2)*x.2+feval(f,A(i,3)*(y+x); c=feval(f,A(i,3)*x; h=-2*c/(b+sign(b)*sqrt(b.2-4*a*c); A(i,4)=A(i,3)+h*(A(i,3)-A(i,2); if(h*(A(i,3)-A(i,2)ep) break end A(i+1,1)=A(i,2); A(i+1,2)=A(i,3); A(i+1,3)=A(i,4);end执行:x0=3;x1=1; x2=1.8;ep=10.-5;f=inline(2(-x)+exp(x)+2*cos(x)-6);muller(f,x0,x1,x2,ep)结果:ans = 3.0000 1.0000 1.8000 1.8218 1.0000 1.8000 1.8218 1.8294 1.8000 1.8218 1.8294 1.82942、用Gauss列主元消去法求解下列方程组 function X=gauss(A,B)%A为系数矩阵%B为值矩阵C=A B;n=length(B);for i=1:n y,c=max(C(i:n,i); D=C(i,:); C(i,:)=C(c+i-1,:); C(c+i-1,:)=D; if C(i,i)=0 break else C(i+1:n,i)=C(i+1:n,i)/C(i,i); C(i+1:n,i+1:n+1)=C(i+1:n,i+1:n+1)-C(i+1:n,i)*C(i,i+1:n+1); endendC(n,n+1)=C(n,n+1)/C(n,n);for i=n-1:-1:1 C(1:i,n+1)=C(1:i,n+1)-C(1:i,i+1)*C(i+1,n+1); C(i,n+1)=C(i,n+1)/C(i,i);endX=C(1:n,n+1);执行:A=4 0 6 0 2;0 1 3 0 2;6 3 19 2 6;0 0 2 5 -5;2 2 6 -5 16;B=12 6 36 2 21;gauss(A,B)结果:ans = 1.0000 1.0000 1.0000 1.0000 1.00003、已知函数表x00.10.20.30.4f（x）0.50.53890.57930.61790.7554分别利用Lagrange和Newton插值法计算的近似值。Lagrange插值法:function yi=lagrange(x,y,xi)%x,y分别为插值节点和节点处的值%xi为被估计函数自变量m=length(x);n=length(y);if(m=n) error(input is wrong);endl=ones(1,n);for i=1:n for j=1:n if i=j continue; else l(1,j)=l(1,j)*(xi-x(1,i)/(x(1,j)-x(1,i); end endend yi=sum(y.*l); 执行:x=0 0.1 0.2 0.3 0.4;y=0.5 0.5389 0.5793 0.6179 0.7554;xi=0.13;lagrange(x,y,xi)结果:ans =0.5530Newton法:function yi=newton(x,y,xi)%x,y分别为插值节点和节点处的值%xi为被估计函数自变量n=length(x);m=length(y);if(m=n) error(input is wrong);endA=;A(1,1:n)=y;for i=1:n-1 for k=1:n-i A(i+1,k)=(A(i,k+1)-A(i,k)/(x(k+i)-x(k); endendB=ones(5,1);for j=2:5 B(j:n,1)=B(j:n,1)*(xi-x(j-1);endyi=sum(B.*A(1:n,1);执行:x=0 0.1 0.2 0.3 0.4;y=0.5 0.5389 0.5793 0.6179 0.7554;xi=0.13;newton(x,y,xi)结果:ans = 0.55304、已知函数及定义区间，将定义区间分成 10等分。利用方法编制三次样条插值函数，满足第一类边界条件：，并输出的值。（注：要求程序具有通用性）function yi=cubicspline(x,y,ydot,xi)%x是向量，是全部的插值点%y是向量，是插值节点出的函数值值%ydot是向量，端点出的导数值%xi为标量%yi为xi出的估计值n=length(x);m=length(y);if n=m error(input is error ,the length of x and y must be equal!); return;endh=zeros(1,n);lamda=ones(1,n);mu=ones(1,n);d=zeros(n,1);for k=1:n-1 h(k)=x(k+1)-x(k); endfor k=1:n-2 lamda(k)=h(k+1)/(h(k)+h(k+1); mu(k)=h(k)/(h(k)+h(k+1); d(k+1)=6/(h(k)+h(k+1)*(y(k+2)-y(k+1)/h(k+1)-(y(k+1)-y(k)/h(k);endA=diag(2*ones(1,n);for i=2:n-1 A(i,i+1)=lamda(i-1); A(i,i-1)=mu(i-1);endA(1,2)=1;A(n,n-1)=1;d(1)=6/h(1)*(y(2)-y(1)/h(1)-ydot(1);d(n)=6/h(n-1)*(ydot(2)-(y(n)-y(n-1)/h(n-1);M=Ad;c=zeros(1,n-1);e=zeros(1,n-1);for i=1:n-1 c(i)=(y(i+1)-y(i)/h(i)-h(i)*(M(i+1)-M(i)/6; e(i)=(y(i)/h(i)-h(i)*M(i)/6)*x(i+1)+(-y(i+1)/h(i)+h(i)*M(i+1)/6)*x(i);endfor k=2:n if x(k-1)=xi&xi=x(k)yi=M(k-1)*