數(shù)值分析編程總結(jié)

1、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); 這個(gè)-1說明行位置-1c2=12*ones(1,n);a2=diag(c2); c3=ones(1,n-1);a3=diag(c3,1); a=a1+a2+a

2、3;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要

3、擬合的函數(shù)為：y=a+bx-cxy 注意不是多項(xiàng)式2. 編程函數(shù)為：function z = erchen(x,y)x1=ones(5,1); A=x1,x,-x.*y; 注意點(diǎn)乘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要擬合的函數(shù)為：

4、y=a+bx+cx2 是多項(xiàng)式2. 編程函數(shù)為：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最小二乘多項(xiàng)式擬合的簡單函數(shù)方法： 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) 要擬

5、合成4次，則2改成4就可以了P =0.8437 0.8642 1.0051 注意此內(nèi)置函數(shù)輸出的結(jié)果c,b,a是反的5 復(fù)合辛普森公式求解積分先定義函數(shù)：function v=f(x) v=sin(x); “若定義有除數(shù)要點(diǎn)除，分母有0時(shí)要特殊定義”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.不動(dòng)點(diǎn)迭代思路

6、不動(dòng)點(diǎn)迭代常常有好幾個(gè)迭代的不動(dòng)點(diǎn)函數(shù),所以要分別定義這些函數(shù)是很困難的,如是乎使用SWITCH內(nèi)置函數(shù)進(jìn)行切換,叫切換函數(shù).1.先定義函數(shù)后進(jìn)行編程的方法先需要定義不動(dòng)點(diǎn)函數(shù):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此函數(shù)的調(diào)用： fixpnt1(f,2,100,1e-5)ans = 13.利用切換函數(shù)SWITCH的方法（多個(gè)不動(dòng)點(diǎn)迭代函數(shù)）function x,it=fix

7、pnt(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

8、,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) br

9、eak;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.先定義函數(shù)后進(jìn)行編程的方法先需要定

10、義不動(dòng)點(diǎn)函數(shù)需要計(jì)算的函數(shù)ffunction v = f( x )v=x5-3*x-10;end需要計(jì)算的函數(shù)的導(dǎo)數(shù)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.先定義函數(shù)后進(jìn)行編程的方法先需要定義不動(dòng)點(diǎn)函數(shù)需要計(jì)算的函數(shù)ffunction v = f(x)v=x2+sin(10*x)-1end需要計(jì)算的函數(shù)的導(dǎo)數(shù)gfunct

11、ion 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:

12、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.先定義函數(shù)后進(jìn)行編程的方法先需要定義函數(shù)需要計(jì)算的函數(shù)ffunction v = f( x )v=x5-3*x-10;end2. 再定義編程:function v

13、 = 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=

14、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;

15、 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. 結(jié)合原點(diǎn)平移反冪法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;whi

16、le (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. 先定義導(dǎo)數(shù)函數(shù)f：function v = f( x,y)v=-y+x+1;end2. 定義歐拉公式編程：function x,y = oula(

17、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. 改進(jìn)的歐拉公式初值問題1. 先定義導(dǎo)數(shù)函數(shù)f：function v = f( x,y)v=-y+x+1;end要先運(yùn)行一下。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(

18、i+1)=0.5*(yc+yp);endend19. 梯形公式初值問題1. 先定義導(dǎo)數(shù)函數(shù)f：function v = f( x,y)v=-y+x+1;end要先運(yùn)行一下。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. 標(biāo)準(zhǔn)四階四段龍格庫塔公式初值問題1. 先定義導(dǎo)數(shù)函數(shù)f：function v = f( x,y)v=y+x;end要先運(yùn)行一下。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*k