数值分析代码

1、第二章 求解线性方程组的数值方法2.1 Gauss消去法function x=gauss(a,b)n=length(b);for k=1:n for i=k+1:n for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a(k,j)/a(k,k); end b(i)=b(i)-a(i,k)*b(k)/a(k,k); endfor i=n:-1:1 %以下为回代过程 s=b(i); if i=n; x(n)=b(n)/a(n,n); else end for j=i+1:n; s=s-a(i,j)*x(j); end x(i)=s/a(i,i);endendend例子：输入矩阵a=2

2、,-1,3;4,2,5;1,2,0; b=1;4;7; 调用函数x=gauss(a,b) 答案为9 -1 -6列主元Gauss消去法function x=gausslie(a,b)n=length(b);for i=1:n max=abs(a(i,i); m=i; for j=i+1:n if maxepsin x0=x1; x1=B*x0+f;endx=x1;例子 A=8,-3,2;4,11,-1;6,3,12 b=20;33;36 x=jacobi (A,b) 答案为3 2 1 Gauss-Seidel迭代法function x=gaussseidel(A,b,p,eps) %p为x的零矩

3、阵,eps为精度L=-tril(A,-1); D=diag(diag(A); U=-triu(A,1); G=(inv(D-L)*U; f=(inv(D-L)*b; x=G*p+f;n=1;while abs(x-p)=eps p=x; x=G*p+f; n=n+1; endn %n为迭代次数的记录例子 A=8,-3,2;4,11,-1;6,3,12 b=20;33;36 p=0;0;0 x=gaussseidel(A,b,p,0.00001) 答案与上一样第三章 非线性方程组的数值解法3.1求非线性方程实根的对分法二分法求根clear all;f=(x) x3-exp(-x); %这是匿名函

4、数，是函数句柄，意思是这个函数已经被写进程序中a=0; b=1; %这个函数是f=x3-e(-x)epsin=1e-8;if f(a)*f(b)epsinif f(a)*f(x1)=epsin x0=x1; x1=1/2*(log(x0)+7); disp(x1); x=max(x1); k=k+1; if k=100 break disp(no solution); endend disp(x);3.3单个非线性方程的Newton法Newton法 书上例子P101，x-e(x)-2=0epsin=1e-5x0=0;x1=x0-(x0+exp(x0)-2)/(1+exp(x0);k=0;whi

5、le abs(x1-x0)=epsin x0=x1 x1=x0-(x0+exp(x0)-2)/(1+exp(x0); k=k+1; if k=100 break disp(no solution); end end disp(x1);Newton正割法（应该就是割线法）clear all;f=(x) x+exp(x)-2;epsin=1e-5;x0=0;x1=1;x2=x1-(f(x1)*(x1-x0)/(f(x1)-f(x0);k=1;while abs(x2-x1)epsin x0=x1; x1=x2; x2=x1-(f(x1)*(x1-x0)/(f(x1)-f(x0); k=k+1;en

6、ddisp(x2);disp(k);例子与上一样第六章插值法6.2拉格朗日插值function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end例子书上P183, x0=0.46 0.47；y0=0.4847 0.4937；x=0.463；y=lagrange(x0,y0,x)6.3Newton插值function

7、 yi=newton(x,y,xi)n=length(x);m=length(y);Y=zeros(n);Y(:,1)=y;for k=1:n-1 for i=1:n-k Y(i,k+1)=(Y(i+1,k)-Y(i,k)/(x(i+k)-x(i); endendyi=0;for i=1:n z=1; for k=1:i-1 z=z*(xi-x(k); end yi=yi+Y(1,i)*z;end例子同上x=0.46 0.47；y=0.4847 0.4937；xi=0.463；yi=newton(x,y,xi)结果一样第八章数值积分复化梯形求积公式function f=fhtx(n,a,b) %f为函数，n为等分数，a,b为端点f=(x) exp(x);h=(b-a)/n;f1=0;f2=f(a)+f(b);for i=0:n-1 x=a+i*h; f1=f1+f(x);endf=(h/2)*(2*f1+f2);复化simpson求积公式function f=simpson(n,a,b);f=(x) exp(x);h=(b-a)/(2*n);f1=0;f2=0;for i=1:n-1 x=a+2*i*h; f1=