数值计算方法matlab程序

1、二分法function x0,k=bisect1(fun1,a,b,ep)if nargin<4ep=1e-5;endfa=feval(fun1,a);fb=feval(fun1,b);if fa*fb>0x0=fa,fb;k=0;return;endk=1;while abs(b-a)/2>epx=(a+b)/2;fx=feval(fun1,x);if fx*fa<0b=x;fb=fx;elsea=x;fa=fx;k=k+1;endendx0=(a+b)/2;>> fun1=inline('x3-x-1');>> x0,k=bi

2、sect1(fun1,1.3,1.4,1e-4)x0 =1.3247k =7>>简单迭代法function x0,k=iterate1(fun1,x0,ep,N)if nargin<4N=500;endif nargin<3ep=1e-5;endx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & k<Nx0=x;x=feval(fun1,x0);k=k+1;endx0=x;if k=Nwarning('已达最大迭代次数')end>> fun1=inline('(x+1)(1/3)'

3、;);>> x0,k=iterate1(fun1,1.5)x0 =1.3247k =7>> fun1=inline('x3-1');>> x0,k=iterate1(fun1,1.5)x0 =Infk =9>>Steffesen加速迭代（简单迭代法的加速）function x0,k=steffesen1(fun1,x0,ep,N)if nargin<4N=500;endif nargin<3ep=1e-5;endx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & k<Nx

4、0=x;y=feval(fun1,x0);z=feval(fun1,y);x=x0-(y-x0)2/(z-2*y+x0);k=k+1;endx0=x;if k=Nwarning('已达最大迭代次数')end>> fun1=inline('(x+1)(1/3)');>> x0,k=steffesen1(fun1,1.5)x0 =1.3247k =3>> fun1=inline('x3-1');>> x0,k=steffesen1(fun1,1.5)x0 =1.3247k =6Newton迭代funct

5、ion x0,k=Newton7(fname,dfname,x0,ep,N)if nargin<5N=500;endif nargin<4ep=1e-5;endx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & k<Nx0=x;x=x0-feval(fname,x0)/feval(dfname,x0);k=k+1;endx0=x;if k=Nwarning('已达最大迭代次数')end>> fname=inline('x-cos(x)');>> dfname=inline(

6、9;1+sin(x)');>> x0,k=Newton7(fname,dfname,pi/4,1e-8)x0 =0.7391k =4非线性方程求根的Matlab函数调用举例：1.求多项式的根： 求f(x)=x3-x-1=0的根：>> roots(1 0 -1 -1)ans =1.3247-0.6624 + 0.5623i-0.6624 - 0.5623i2.求一般函数的根>> fun=inline('x*sin(x2-x-1)','x')fun = Inline function: fun(x) = x*sin(x2-

7、x-1)>> fplot(fun,-2 0.1);grid on>> x=fzero(fun,-2,-1)x = -1.5956>> x=fzero(fun,-1 -0.1)x = -0.6180x,f,h=fsolve(fun,-1.6)x = -1.5956f = 1.4909e-009h = 1（h>0表示收敛，h<0表示发散，h=0表示已达到设定的计算函数值的最大次数）第三章：线性方程组的数值解法1. 高斯消元法function A,x=gauss3(A,b)%本算法用顺序高斯消元法求解线性方程组n=length(b);A=A,b;for

8、 k=1:n-1 A(k+1):n,(k+1):(n+1)=A(k+1):n,(k+1):(n+1)-A(k+1):n,k)/A(k,k)*A(k,(k+1):(n+1); A(k+1):n,k)=zeros(n-k,1); A;endx=zeros(n,1); %上面为消元过程 x(n)=A(n,n+1)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,n+1)-A(k,(k+1):n)*x(k+1:n)/A(k,k); end %上面为回代过程 >> A=2 3 4;3 5 2;4 3 30;>> b=6,5,32'b = 6 5 32&g

9、t;> A,x=gauss3(A,b)A = 2.0000 3.0000 4.0000 6.0000 0 0.5000 -4.0000 -4.0000 0 0 -2.0000 -4.0000x = -13 8 2列选主元的高斯消元法：function A,x=gauss5(A,b)%本算法用列选主元的高斯消元法求解线性方程组n=length(b);A=A,b;for k=1:n-1 %选主元 ap,p=max(abs(A(k:n,k); p=p+k-1; if p>k t=A(k,:); A(k,:)=A(p,:); A(p,:)=t; end %消元 A(k+1):n,(k+1)

10、:(n+1)=A(k+1):n,(k+1):(n+1)-A(k+1):n,k)/A(k,k)*A(k,(k+1):(n+1); A(k+1):n,k)=zeros(n-k,1); end%回代 x=zeros(n,1); x(n)=A(n,n+1)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,n+1)-A(k,(k+1):n)*x(sk+1:n)/A(k,k);end>> A=2 3 4;3 5 2;4 3 30; b=6,5,32'>> A,x=gauss5(A,b)A = 4.0000 3.0000 30.0000 32.0000 0

11、2.7500 -20.5000 -19.0000 0 0 0.1818 0.3636x = -13 8 2三角分解法：Doolittle 分解function L,U=doolittle1(A)n=length(A);U=zeros(n);L=eye(n);U(1,:)=A(1,:);L(2:n,1)=A(2:n,1)/U(1,1);for k=2:n U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n); L(k+1:n,k)=A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,n)/U(k,k);Endy=zeros(n,1);x=y;y(1)

12、=b(1);for i=2:n y(i)=b(i)-L(i,1:i-1)*y(1:i-1);endx(n)=y(n)/U(n,n);for i=n-1:-1:1 x(i)=(y(i)-U(i,i+1:n)*x(i+1:n)/U(i,i);end>> A=1 2 3;2 5 2 ;3 1 5;b=14 18 20'>> L,U,x=doolittle1(A,b)L = 1 0 0 2 1 0 3 -8 1U = 1 2 3 0 1 -4 0 0 -36x = 2.8333 1.33332.8333 平方根法：function L,x=choesky3(A,b)n=

13、length(A);L=zeros(n);L(:,1)=A(:,1)/sqrt(A(1,1);for k=2:n L(k,k)=A(k,k)-L(k,1:k-1)*L(k,1:k-1)' L(k,k)=sqrt(L(k,k); for i=k+1:n L(i,k)=(A(i,k)-L(i,1:k-1)*L(k,1:k-1)')/L(k,k); endend y=zeros(n,1);x=y;y(1)=b(1)/L(1,1);for i=2:n y(i)=(b(i)-L(i,1:i-1)*y(1:i-1)/L(i,i);endx(n)=y(n)/L(n,n);for i=n-1:

14、-1:1 x(i)=(y(i)-L(i+1:n,i)'*x(i+1:n)/L(i,i);end>> A=4 -1 1;-1 4.25 2.75;1 2.75 3.5A = 4.0000 -1.0000 1.0000 -1.0000 4.2500 2.7500 1.0000 2.7500 3.5000>> b=4 6 7.25'b = 4.0000 6.0000 7.2500L,x=choesky3(A,b)L = 2.0000 0 0 -0.5000 2.0000 0 0.5000 1.5000 1.0000x = 1 1 1>>迭代法求方程

15、组的解Jacobi迭代法：function x,k=jacobi2(a,b,x0,ep,N)%本算法用Jacobi迭代求解ax=b,用分量形式n=length(b);k=0;if nargin<5 N=500;endif nargin<4 ep=1e-5;endif nargin<3 x0=zeros(n,1); y=zeros(n,1);endx=x0;x0=x+2*ep;while norm(x-x0,inf)>ep & k<N k=k+1; x0=x; for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-a

16、(i,j)*x0(j); end end if abs(a(i,i)<1e-10|k=N warning('a(i,i) is too small'); return end y(i)=y(i)/a(i,i); end x=y; enda=4 3 0;3 4 -1; 0 -1 4;b=24 30 -24'；x,k=jacobi2(a,b)x = 3.0000 4.0000 -5.0000k =59Gauss-seidel迭代法：function x,k=gaussseide2(a,b,x0,ep,N)%本算法用Gauss-seidel迭代求解ax=b,用分量形式n

17、=length(b);k=0;if nargin<5 N=500;endif nargin<4 ep=1e-5;endif nargin<3 x0=zeros(n,1); y=zeros(n,1);endx=x0;x0=x+2*ep;while norm(x-x0,inf)>ep & k<N k=k+1; x0=x; y=x; for i=1:n z(i)=b(i); for j=1:n if j=i z(i)=z(i)-a(i,j)*x(j); end end if abs(a(i,i)<1e-10|k=N warning('a(i,i)

18、is too small'); return end z(i)=z(i)/a(i,i); x(i)=z(i); end endx,k=gaussseide2(a,b)x = 3.0000 4.0000 -5.0000k = 25最速下降法function x,k=zuisuxiajiang(A,b,x0,ep,N)%本算法用最速下降算法求解正定方程组Ax=b,n=length(b);if nargin<5 N=500;endif nargin<4 ep=1e-8;endif nargin<3 x0=ones(n,1); endx=x0;x0=x+2*ep;r=b-A*

19、x;d=r;k=0;while norm(x-x0,inf)>ep & k<N k=k+1; x0=x;lamda=(d'*d)/(d'*A*d);x=x0+lamda*d;r=b-A*x;d=r;endif k=N warning('已达最大迭代次数')end共轭梯度算法function x,k=gongertidufa(A,b,x0,ep,N)%本算法用共轭梯度算法求解正定方程组Ax=b,n=length(b);if nargin<5 N=500;endif nargin<4 ep=1e-8;endif nargin<3

20、 x0=ones(n,1); endx=x0;x0=x+2*ep;r=b-A*x;d=r;k=0;while norm(x-x0,inf)>ep & k<N k=k+1; x0=x;lamda=(r'*r)/(d'*A*d);r1=r;x=x0+lamda*d;r=b-A*x;beta=(r'*r)/(r1'*r1);d=r+beta*d;endif k=N warning('已达最大迭代次数')end 常微分方程数值解function x,y=Euler1(fun,xspan,y0,h)%本算法用欧拉格式计算微分方程y

21、9;=f(x,y)的解。 %fun为右端函数，xspan为求解的自变量区间，yo为初值，h为步长。% 输出向量x以及对应的函数值向量y.x=xspan(1):h:xspan(2);y(1)=y0;n=length(x);for k=1:n-1 y(k+1)=y(k)+h*feval(fun,x(k),y(k);endx=x'y=y'fun=inline('y-2*x/y');x,y1=Euler1(fun,0,1,1,0.1)x = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1y1 = 1 1.1 1.1918 1.2774 1

22、.3582 1.4351 1.509 1.5803 1.6498 1.7178 1.7848function x,y=Euler2(fun,xspan,y0,h)%本算法用改进的欧拉格式计算微分方程y'=f(x,y)的解。 %fun为右端函数，xspan为求解的自变量区间，yo为初值，h为步长。% 输出向量x以及对应的函数值向量y.x=xspan(1):h:xspan(2);y(1)=y0;n=length(x);for k=1:n-1 k1=feval(fun,x(k),y(k); y(k+1)=y(k)+h*k1; k2=feval(fun,x(k+1),y(k+1); y(k+1)=y(k)+h*(k1+k2)/2;endx=x'y=y'x,y2=Euler2(fun,0,1,1,0.1)x = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1y2 = 1 1.0959 1.1841 1.2662 1.3434 1.4164 1.4