数值计算与符号计算.pptVIP

再调用dblquad函数求解。 global ki;ki=0; I=dblquad(fxy, -2, 2, -1, 1) ki 程序运行结果如下： I= 1.57449318974494 ki= 1038 4.3 线性方程组求解 对于线性方程组Ax=b，可以利用左除运算符“\”求解： x=A\b 即x=inv(A)*b 如果矩阵A是奇异的或接近奇异的，则MATLAB会给出警告信息。 MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式如下。 [L, U]=lu(X)，产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。 [L, U, P]=lu(X)，产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。 实现LU分解后，线性方程组Ax=b的解可表示为x=U\(L\b)或x=U\(L\Pb)，这样可以大大提高运算速度。 MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式如下。 [Q, R]=qr(X)，产生一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。 [Q, R, E]=qr(X)，产生一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。 实现QR分解后，线性方程组Ax=b的解可表示为x=R\(Q\b)或x=E(R\(Q\b))。 MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式如下。 R=chol(X)，产生一个上三角阵R，使RR=X。若X为非对称正定的，则输出一个出错信息。 [R, p]=chol(X)，这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足RR=X(1:q, 1:q)。 实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x=R\(R\b)。 例4.12 用直接法求解下列线性方程组。 MATLAB程序如下： A=[1/2 1/3 1/4;1/3 1/4 1/5;1/4 1/5 1/6] B=[0.95 0.67 0.52] X=A\B或X=inv(A)*B 程序运行结果如下： X= 1.2000 0.6000 0.6000 例4.13 分别用LU分解、QR分解和Cholesky分解求解例4.12中的线性方程组。 MATLAB程序如下： A=[1/2 1/3 1/4;1/3 1/4 1/5;1/4 1/5 1/6] B=[0.95 0.67 0.52] [L, U]=lu(A); X=U\(L\B) %LU分解求解 [Q, R]=qr(A); X=R\(Q\B) %QR分解求解 R=chol(A) X=R\(R\B) %Cholesky分解求解 所得结果均为 X= 1.2000 0.6000 0.6000 Jacobi迭代法的MATLAB函数文件Jacobi.m如下： function [y, n]=jacobi(A, b, x0, eps) if nargin==3 eps=1.0e-6; elseif nargin3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A, -1); %求A的下三角阵 U=-triu(A, 1); %求A的上三角阵 B=D\(L+U); f=D\b; y=B*x0+f; n=1; %迭代次数 while norm(y-x0)=eps x0=y; y=B*x0+f; n=n+1; end 例4.14 用Jacobi迭代法求解线性方程组。设迭代初值为0，迭代精度为10-6。 在命令中调用函数文件Jacobi.m，命令如下： A=[9, -1, 1;-1, 10, -2; -2, 1, 10]; b=[10, 7, 6]; [x, n]=jacobi(A, b, [0, 0, 0], 1.0e-6) MATLAB返回 x= 1.1365 0.9599 0.7313 n= 10 输入x=inv(A)*b，所得结果相同。 4.4 常微分方程的数值求解 MATLAB提供了求常微分方程数值解的函数，一般调用格式如下： [t, y]=ode23(fname, tspan, y0), 二阶、三阶龙格-库塔法 [t, y]=ode45(fname, tspan, y0)，四阶、五阶龙格-库塔法 其中，fn