数值分析编程总结

个人精品文档资料1. LU分解：L U=lu(A);2. 追赶法function x=zhuiganfa(A,b)n,n=size(A); for i=1:n if(i=1) l(i)=A(i,i); y(i)=b(i)/l(i); else l(i)=A(i,i)-A(i,i-1)*u(i-1); y(i)=(b(i)-y(i-1)*A(i,i-1)/l(i); end if(in=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. 拉格朗日插值function yh=lage(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);for i=1:n xp=x(1:i-1 i+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 最小二乘法1. 最小二乘的xi和yi为：xi19531964198219902000yi5.826.9510.0811.3412.66要拟合的函数为：y=a+bx-cxy 注意不是多项式2. 编程函数为：function z = erchen(x,y)x1=ones(5,1); A=x1,x,-x.*y; 注意点乘z=Ay; 注意左除a=z(1);b=z(2);c=z(3);end 输入： x=1953 1964 1982 1990 2000; y=5.82 6.95 10.08 11.34 12.66; erchen(x,y)ans = 2.9456 =a -0.0014 =b -0.0005 =c1. 最小二乘的xi和yi为：xi00.250.50.751yi11.2841.64872.11702.7183要拟合的函数为：y=a+bx+cx2 是多项式2. 编程函数为：function z = erchen2(x,y)x1=ones(5,1); A=x1,x,x.2;z=Ay;a=z(1);b=z(2);c=z(3);end输入： x=0 0.25 0.5 0.75 1.00; y=1.00 1.284 1.6487 2.1170 2.7183; erchen2(x,y)ans = 1.0051 0.8642 0.8437最小二乘多项式拟合的简单函数方法： x=0 0.25 0.5 0.75 1.00; y=1.00 1.284 1.6487 2.1170 2.7183; P=polyfit(x,y,2) 要拟合成4次，则2改成4就可以了P =0.8437 0.8642 1.0051 注意此内置函数输出的结果c,b,a是反的5 复合辛普森公式求解积分先定义函数：function v=f(x) v=sin(x); “若定义有除数要点除，分母有0时要特殊定义”end定义程序：function I=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.先定义函数后进行编程的方法先需要定义不动点函数:function v=f(x)v=x3-x-1;end再定义编程:function it,x=fixpnt1(f,a,maxit,tol)it=0;x=feval(f,a);while ittol, 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)switch np, case 1, phi=inline(3*x+10)(1/5); case 2, phi=inline(sin(10*x)+2*cos(x)-3); case 3, phi=inline(3-atan(x); case 4, phi=inline(-2-1/log(x2+x+1);endit=0;x=phi(a);while ittol, it=it+1; a=x; x=phi(a); end使用与输入： fixpnt(2,1,100,1e-5)ans = -4.26967. 雅可比迭代function x it=jacobi(A,b,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);x=zeros(size(b);for it=1:500 x=D(b+L*x+U*x);error=norm(b-A*x)/norm(b);if (errortol) break;endend8. 高斯迭代function x it=gaosi(A,b,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);x=zeros(size(b);for it=1:500 x=(D-L)(b+U*x);error=norm(b-A*x)/norm(b);if (errortol) break;endend9. SOR迭代function x it=SOR(A,b,w,tol)D=diag(diag(A);L=D-tril(A);U=D-triu(A);x=zeros(size(b);for it=1:500 x=(D-w*L)(w*b+(1-w)*D*x+w*U*x);error=norm(b-A*x)/norm(b);if (errortol, it=it+1; x=a/2+b/2; fx=feval(f,x);if sign(fx)=sign(fa), a=x;fa=fx;else b=x;fb=fx;end end11. 牛顿法：1.先定义函数后进行编程的方法先需要定义不动点函数需要计算的函数ffunction v = f( x )v=x5-3*x-10;end需要计算的函数的导数gfunction v = g( x )v=5*x4-3;end2. 再定义编程:function v = newton(a,f,g,maxit,tol)it=0;x=a;while ittol, it=it+1; x=x-feval(f,x)/feval(g,x);endv=x,it;end12. 牛顿下山法：1.先定义函数后进行编程的方法先需要定义不动点函数需要计算的函数ffunction v = f(x)v=x2+sin(10*x)-1end需要计算的函数的导数gfunction v =g(x)v = 2*x+10*cos(10*x)end2. 再定义编程1:function v = newtonxiashang(x0,f,g,maxit,tol)x=x0;it=0;while ittol, it=it+1; d=-feval(f,x)/feval(g,x); lambda=1; isdone=0; while isdone, xn=x+lambda*d; if abs(feval(f,xn)abs(feval(f,x), isdone=1; else lambda=lambda*0.5; end end x=xn;endv=x,it;end3. 再定义编程2:function v = newtonxiashang2(x0,f,g,maxit,tol)x=x0;it=0;while ittol, it=it+1; d=-feval(f,x)/feval(g,x); lambda=1; while abs(feval(f,x+lambda*d)=abs(feval(f,x), lambda=0.5*lambda; end x=x+lambda*d;end v=x,it;end13. 割线法1.先定义函数后进行编程的方法先需要定义函数需要计算的函数ffunction v = f( x )v=x5-3*x-10;end2. 再定义编程:function v = gexian(a,b,f,maxit,tol)it=0;x0=a;x=b;while ittol, 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.7226 7.000014. 乘幂法function t,y = chenmifa( a,xinit,ep )v0=xinit;tv ti= max(abs(v0);lam0=v0(ti);u0=v0/lam0;flag=0;while (flag=0) v1=a*u0; tv ti= max(abs(v1); lam1=v1(ti); u0=v1/lam1; if (abs(lam0-lam1) a=2,3,2;10,3,4;3,6,1; xinit=1 1 1; ep=0.0001; t y=chenmifa(a,xinit,ep)t = 11.0000y = 0.5000 1.0000 0.750015. 反幂法function t y =fanmifa( a,xinit,ep )v0=xinit;tv ti= max(abs(v0);lam0=v0(ti);u0=v0/lam0;flag=0;while (flag=0) v1=a-1*u0; tv ti= max(abs(v1); lam1=v1(ti); u0=v1/lam1; if (abs(lam0-1-lam1-1) a=12 6 -6;6 16 2;-6 2 16;xinit=1 -0.5 0.5； t y =fanmifa( a,xinit,0.0001 )t = 4.4560y = 1.0000 -0.62870.628716. 结合原点平移反幂法function t y =yuandian( a,p,xinit,ep )n n=size(a);v0=xinit;tv ti= max(abs(v0);lam0=v0(ti);u0=v0/lam0;flag=0;while (flag=0) v1=(a-p*eye(n)-1*u0; tv ti= max(abs(v1); lam1=v1(ti); u0=v1/lam1; if (abs(lam0-1-lam1-1) a=6 2 1;2 3 1;1 1 1; p=6; ep=0.0001; xinit=1 1 1; t y=yuandian(a,p,xinit,ep)t = 7.2880y = 1.0000 0.5229 0.242217. 欧拉公式初值问题1. 先定义导数函数f：function v = 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;for i=1:n; y(i+1)=y(i)+h*feval(f,x(i),y(i);endend18. 改进的欧拉公式初值问题1. 先定义导数函数f：function v = 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;for i=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：function v = f( x,y)v=-y+x+1;end要先运行一下。2. 定义欧拉公式编程：function x,y = tixing( f,y0,a,b,n )y(1)=y0;h=(b-a)/n;x=a:h:b;for i=1:n; y(i+1)=1/(0.5*h+1)*(y(i)+0.5*h*(feval(f,x(i),y(i)+x(i+1)+1);endend20. 标准四阶四段龙格库塔公式初值问题1. 先定义导数函数f：function v = f( x,y)v=y+x;end要先运行一下。2. 定义欧拉公式编程：function x,y = longgekuta(f,y0,a,b,n)y(1)=y0;h=(b-a)/n;x=a:h:b;for i=1:n, k1=h*feval(f,x(i),y(i); k2=h*feval(f,x(i)+0.5*h,y(i)+0.5*k1); k3=h*feval(f,x(i)+0.5*h,y(i)+0.5*k2); k4=h*feval(f,x(i)+h,y(i)+k3); y(i+1)=y(i)+(k1+2*k2+2*k3+