数值分析编程总结

1.LU分解：[LU]=lu(A);2.追赶法functionx=zhuiganfa(A,b)[n,n]=size(A);fori=1:nif(i==1)l(i)=A(i,i);y(i)=b(i)/l(i);elsel(i)=A(i,i)-A(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);endif(i<n)u(i)=A(i,i+1)/l(i);endendx(n)=y(n)forj=n-1:-1:1x(j)=y(j)-u(j)*x(j+1);end数值试验：n=101;a=1210………………..….01121………………..…001121………..…0001121……01112关键是如何定义上述矩阵：>>n=101;c1=ones(1,n-1);a1=diag(c1,-1);这个-1说明行位置-1c2=12*ones(1,n);a2=diag(c2);c3=ones(1,n-1);a3=diag(c3,1);a=a1+a2+a3;3.拉格朗日插值functionyh=lage(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);fori=1:nxp=x([1:i-1i+1:n]);yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));endend>>x=[11,12];>>y=[2,4];>>xh=[11.75];>>lage(x,y,xh)ans=3.50004最小二乘法最小二乘的xi和yi为：xi19531964198219902000yi5.826.9510.0811.3412.66要拟合的函数为：y=a+bx-cxy注意不是多项式编程函数为：functionz=erchen(x,y)x1=ones(5,1);A=[x1,x,-x.*y];注意点乘z=A\y;注意左除a=z(1);b=z(2);c=z(3);end输入：≻≻x=[19531964198219902000]';≻≻y=[5.826.9510.0811.3412.66]';≻≻erchen(x,y)ans=2.9456=a-0.0014=b-0.0005=c最小二乘的xi和yi为：xi00.250.50.751yi11.2841.64872.11702.7183要拟合的函数为：y=a+bx+cx2是多项式编程函数为：functionz=erchen2(x,y)x1=ones(5,1);A=[x1,x,x.^2];z=A\y;a=z(1);b=z(2);c=z(3);end输入：≻≻x=[00.250.50.751.00]';≻≻y=[1.001.2841.64872.11702.7183]';≻≻erchen2(x,y)ans=1.00510.86420.8437最小二乘多项式拟合的简单函数方法：≻≻x=[00.250.50.751.00]';≻≻y=[1.001.2841.64872.11702.7183]';≻≻P=polyfit(x,y,2)要拟合成4次，则2改成4就可以了P=0.84370.86421.0051注意此内置函数输出的结果c,b,a是反的5复合辛普森公式求解积分先定义函数：functionv=f(x)v=sin(x);“若定义有除数要点除，分母有0时要特殊定义”end定义程序：functionI=fsps(f,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval(f,x);I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));end>>fsps('f',0,1,4)ans=0.45976.不动点迭代思路不动点迭代常常有好几个迭代的不动点函数,所以要分别定义这些函数是很困难的,如是乎使用SWITCH内置函数进行切换,叫切换函数.1.先定义函数后进行编程的方法先需要定义不动点函数:functionv=f(x)v=x^3-x-1;end再定义编程:function[it,x]=fixpnt1(f,a,maxit,tol)it=0;x=feval(f,a);whileit<=maxit&abs(x-a)>tol,it=it+1;a=x;x=feval(f,a);end此函数的调用：>>fixpnt1('f',2,100,1e-5)ans=13.利用切换函数SWITCH的方法（多个不动点迭代函数）function[x,it]=fixpnt(np,a,maxit,tol)switchnp,case1,phi=inline('(3*x+10)^(1/5)');case2,phi=inline('sin(10*x)+2*cos(x)-3');case3,phi=inline('3-atan(x)');case4,phi=inline('-2-1/log(x^2+x+1)');endit=0;x=phi(a);whileit<=maxit&abs(x-a)>tol,it=it+1;a=x;x=phi(a);end使用与输入：>>fixpnt(2,1,100,1e-5)ans=-4.26967.雅可比迭代function[xit]=jacobi(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));forit=1:500x=D\(b+L*x+U*x);error=norm(b-A*x)/norm(b);if(error<tol)break;endend8.高斯迭代function[xit]=gaosi(A,b,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));forit=1:500x=(D-L)\(b+U*x);error=norm(b-A*x)/norm(b);if(error<tol)break;endend9.SOR迭代function[xit]=SOR(A,b,w,tol)D=diag(diag(A));L=D-tril(A);U=D-triu(A);x=zeros(size(b));forit=1:500x=(D-w*L)\(w*b+(1-w)*D*x+w*U*x);error=norm(b-A*x)/norm(b);if(error<tol)break;endend10.二分法：1.先要定义所求的函数:functionv=f(x)v=x^3-x-1;end2.二分法程序如下:function[x,it]=erfengfa(a,b,f,tol)fa=feval(f,a);fb=feval(f,b);it=0;whileabs(b-a)>tol,it=it+1;x=a/2+b/2;fx=feval(f,x);ifsign(fx)==sign(fa),a=x;fa=fx;elseb=x;fb=fx;endend11.牛顿法：1.先定义函数后进行编程的方法先需要定义不动点函数需要计算的函数ffunctionv=f(x)v=x^5-3*x-10;end需要计算的函数的导数gfunctionv=g(x)v=5*x^4-3;end2.再定义编程:functionv=newton(a,f,g,maxit,tol)it=0;x=a;whileit<=maxit&abs(feval(f,x))>tol,it=it+1;x=x-feval(f,x)/feval(g,x);endv=[x,it];end12.牛顿下山法：1.先定义函数后进行编程的方法先需要定义不动点函数需要计算的函数ffunctionv=f(x)v=x^2+sin(10*x)-1end需要计算的函数的导数gfunctionv=g(x)v=2*x+10*cos(10*x)end2.再定义编程1:functionv=newtonxiashang(x0,f,g,maxit,tol)x=x0;it=0;whileit<=maxit&abs(feval(f,x))>tol,it=it+1;d=-feval(f,x)/feval(g,x);lambda=1;isdone=0;while~isdone,xn=x+lambda*d;ifabs(feval(f,xn))<abs(feval(f,x)),isdone=1;elselambda=lambda*0.5;endendx=xn;endv=[x,it];end3.再定义编程2:functionv=newtonxiashang2(x0,f,g,maxit,tol)x=x0;it=0;whileit<=maxit&abs(feval(f,x))>tol,it=it+1;d=-feval(f,x)/feval(g,x);lambda=1;whileabs(feval(f,x+lambda*d))>=abs(feval(f,x)),lambda=0.5*lambda;endx=x+lambda*d;endv=[x,it];end13.割线法1.先定义函数后进行编程的方法先需要定义函数需要计算的函数ffunctionv=f(x)v=x^5-3*x-10;end2.再定义编程:functionv=gexian(a,b,f,maxit,tol)it=0;x0=a;x=b;whileit<=maxit&abs(feval(f,x))>tol,it=it+1;xt=x-feval(f,x)/(feval(f,x)-feval(f,x0))*(x-x0);x0=x;x=xt;endv=[x,it];end输入:>>gexian(1,2,'f',100,1e-5)ans=1.72267.000014.乘幂法function[t,y]=chenmifa(a,xinit,ep)v0=xinit;[tvti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while(flag==0)v1=a*u0;[tvti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if(abs(lam0-lam1)<=ep)flag=1;endlam0=lam1;endt=lam1;y=u0;end运行及其结果:>>a=[2,3,2;10,3,4;3,6,1];>>xinit=[111]';>>ep=0.0001;>>[ty]=chenmifa(a,xinit,ep)t=11.0000y=0.50001.00000.750015.反幂法function[ty]=fanmifa(a,xinit,ep)v0=xinit;[tvti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while(flag==0)v1=a^-1*u0;[tvti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if(abs(lam0^-1-lam1^-1)<=ep)flag=1;endlam0=lam1;endt=lam1^-1;y=u0;end运行提示：>>a=[126-6;6162;-6216];>>xinit=[1-0.50.5]'；>>[ty]=fanmifa(a,xinit,0.0001)t=4.4560y=1.0000-0.62870.628716.结合原点平移反幂法function[ty]=yuandian(a,p,xinit,ep)[nn]=size(a);v0=xinit;[tvti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;flag=0;while(flag==0)v1=(a-p*eye(n))^-1*u0;[tvti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if(abs(lam0^-1-lam1^-1)<=ep)flag=1;endlam0=lam1;endt=lam1^-1+p;y=u0;end运行结果:>>a=[621;231;111];>>p=6;>>ep=0.0001;>>xinit=[111]';>>[ty]=yuandian(a,p,xinit,ep)t=7.2880y=1.00000.52290.242217.欧拉公式初值问题1.先定义导数函数f：functionv=f(x,y)v=-y+x+1;end2.定义欧拉公式编程：function[x,y]=oula(f,y0,a,b,n)y(1)=y0;h=(b-a)/n;x=a:h:b;fori=1:n;y(i+1)=y(i)+h*feval(f,x(i),y(i));endend18.改进的欧拉公式初值问题1.先定义导数函数f：functionv=f(x,y)v=-y+x+1;end要先运行一下。2.定义欧拉公式编程：function[x,y]=gaijingoula(f,y0,a,b,n)y(1)=y0;h=(b-a)/n;x=a:h:b;fori=1:n;yp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=0.5*(yc+yp);endend19.梯形公式初值问题1.先定义导数函数f：functionv=f(x,y)v=-y+x+1;end要先运行一下。2.定义欧拉公式编程：function[x,y]=tixing(