搜索
第二章求解线性方程组的数值方法2.1Gauss消去法functionx=gauss(a,b)n=length(b);fork=1:nfori=k+1:nforj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j)/a(k,k);endb(i)=b(i)-a(i,k)*b(k)/a(k,k);endfori=n:-1:1%以下为回代过程s=b(i);ifi==n;x(n)=b(n)/a(n,n);elseendforj=i+1:n;s=s-a(i,j)*x(j);endx(i)=s/a(i,i);endendend例子：输入矩阵a=[2,-1,3;4,2,5;1,2,0];b=[1;4;7];调用函数x=gauss(a,b)答案为9-1-6列主元Gauss消去法functionx=gausslie(a,b)n=length(b);fori=1:nmax=abs(a(i,i));m=i;forj=i+1:nifmaxabs(a(j,i))max=abs(a(j,i));m=j;endendendif(m~=i)fork=i:nc(k)=a(i,k);a(i,k)=a(m,k);a(m,k)=c(k);endb=b(i);b(i)=b(m);b(m)=b;end%以上是搜索绝对值最大元和交换方程的代码（最大元素法）fork=1:nfori=k+1:nforj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j)/a(k,k);endb(i)=b(i)-a(i,k)*b(k)/a(k,k);endfori=n:-1:1s=b(i);ifi==n;x(n)=b(n)/a(n,n);elseendforj=i+1:n;s=s-a(i,j)*x(j);endx(i)=s/a(i,i);endendend%后面这段就是Gauss消去法的代码例子可以用上面的，结果是一样的2.2迭代法需注意迭代格式是怎样推出来的（逐次逼近法），注意L-U分解是怎样实现的Jacobi迭代法functionx=jacobi(A,b)n=length(A);D=zeros(n,n);L=zeros(n,n);U=zeros(n,n);fori=1:n;D(i,i)=A(i,i);endfori=2:nforj=1:i-1L(i,j)=-A(i,j);U(j,i)=-A(j,i);endendB=inv(D)*(L+U);f=inv(D)*b;x0=ones(n,1);x1=B*x0+f;epsin=1e-5;whilemax(abs(x1-x0))epsinx0=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)答案为321Gauss-Seidel迭代法functionx=gaussseidel(A,b,p,eps)%p为x的零矩阵,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;whileabs(x-p)=epsp=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求非线性方程实根的对分法二分法求根clearall;f=@(x)x^3-exp(-x);%这是匿名函数，@是函数句柄，意思是这个函数已经被写进程序中a=0;b=1;%这个函数是f=x^3-e^(-x)epsin=1e-8;iff(a)*f(b)0x1=(a+b)/2;endwhileabs(x1-b)epsiniff(a)*f(x1)0;b=x1;x1=(a+b)/2;elsea=x1;x1=(a+b)/2;endendxstar=x1;disp(x1);这个不需要调用，直接运行就有结果了3.2单个非线性方程的迭代法注：迭代格式必须是收敛的epsin=1e-5;x0=3.5;x1=1/2*(log(x0)+7);k=0;whileabs(x1-x0)=epsinx0=x1;x1=1/2*(log(x0)+7);disp(x1);x=max(x1);k=k+1;ifk==100breakdisp('nosolution');endenddisp(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;whileabs(x1-x0)=epsinx0=x1x1=x0-(x0+exp(x0)-2)/(1+exp(x0));k=k+1;ifk==100breakdisp('nosolution');endenddisp(x1);Newton正割法（应该就是割线法）clearall;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;whileabs(x2-x1)epsinx0=x1;x1=x2;x2=x1-(f(x1)*(x1-x0))/(f(x1)-f(x0));k=k+1;enddisp(x2);disp(k);例子与上一样第六章插值法6.2拉格朗日插值functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end例子书上P183,x0=[0.460.47]；y0=[0.48470.4937]；x=0.463；y=lagrange(x0,y0,x)6.3Newton插值functionyi=newton(x,y,xi)n=length(x);m=length(y);Y=zeros(n);Y(:,1)=y';fork=1:n-1fori=1:n-kY(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));endendyi=0;fori=1:nz=1;fork=1:i-1z=z*(xi-x(k));endyi=yi+Y(1,i)*z;end例子同上x=[0.460.47]；y=[0.48470.4937]；xi=0.463；yi=newton(x,y,xi)结果一样第八章数值积分复化梯形求积公式functionf=fhtx(n,a,b)%f为函数，n为等分数，a,b为端点f=@(x)exp(x);h=(b-a)/n;f1=0;f2=f(a)+f(b);fori=0:n-1x=a+i*h;f1=f1+f(x);endf=(h/2)*(2*f1+f2);复化simpson求积公式functionf=simpson(n,a,b);f=@(x)exp(x);h=(b-a)/(2*n);f1=0;f2=0;fori=1:n-1x=a+2*i*h;f1=f1+f(x);endfori=1:nx=a+(2*i-1)*h;f2=f2+f(x);endsnf=(h/3)*(f(a)+2*f1+4*f2+f(b))调用的时候需要输入等分数n，端点a,b。函数也是已经写进代码了，可以根据需要更改。