《刘卫国全套配套课件matlab程序设计与应用第3版第8章matlab方程数值求解源程序》由会员分享，可在线阅读，更多相关《刘卫国全套配套课件matlab程序设计与应用第3版第8章matlab方程数值求解源程序（7页珍藏版）》请在金锄头文库上搜索。
1、第8章 MATLAB方程数值求解例8-1 用直接解法求解下列线性方程组。程序如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;x=Ab例8-2 用LU分解求解例8-1中的线性方程组。程序如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;L,U=lu(A);x=U(Lb)例8-3 用QR分解求解例8-1中的线性方程组。程序如下：A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;Q,R=qr(A);x=R(Qb)例8-4 用Chol
2、esky分解求解例8-1中的线性方程组。命令如下： A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; R=chol(A)Jacobi迭代法的MATLAB函数文件jacobi.m如下：function y,n=jacobi(A,b,x0,ep)if nargin=3 ep=1.0e-6;elseif nargin=ep x0=y; y=B*x0+f; n=n+1;end例8-5 用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在程序中调用函数文件jacobi.m，程序如下：A=10,-1,0;-1,10,-2;0
3、,-2,10;b=9,7,6;x,n=jacobi(A,b,0,0,0,1.0e-6)Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：function y,n=gauseidel(A,b,x0,ep)if nargin=3 ep=1.0e-6;elseif nargin=ep x0=y; y=G*x0+f; n=n+1;end例8-6 用Gauss-Serdel迭代法求解例8-5中的线性方程组。在程序中调用函数文件gauseidel.m，程序如下：A=10,-1,0;-1,10,-2;0,-2,10;b=9,7,6;x,n=gauseidel(A,b,0,0,0
4、,1.0e-6)例8-7 分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。命令如下： a=1,2,-2;1,1,1;2,2,1; b=9;7;6; x,n=jacobi(a,b,0;0;0) x,n=gauseidel(a,b,0;0;0)有了上面这些讨论，下面设计一个求解线性方程组的函数文件line_solution.m。function x,y=line_solution(A,b)m,n=size(A);y=;if norm(b)0 %非齐次方程组 if rank(A)=rank(A,b) if rank(A)=n %有唯一解 disp(原方程组有唯一
5、解x) x=Ab; else %方程组有无穷多个解，基础解系 disp(原方程组有无穷个解，特解为x，其齐次方程组的基础解系为y) x=Ab; y=null(A,r); end else disp(方程组无解) %方程组无解 x=; endelse %齐次方程组 disp(原方程组有零解x) x=zeros(n,1); %0解 if rank(A) fzero(fz,-5) %以-5作为迭代初值 fzero(fz,1) %以1作为迭代初值例8-11 求下列方程组在（1，1，1）附近的解并对结果进行验证。首先建立函数文件myfun.m。function F=myfun(X)x=X(1);y=X(
6、2);z=X(3);F(1)=sin(x)+y+z2*exp(x);F(2)=x+y+z;F(3)=x*y*z;在给定的初值x0=1，y0=1，z0=1下，调用fsolve函数求方程的根，命令如下： option=optimset(Display,off); X=fsolve(myfun,1,1,1,option) q=myfun(X)例8-12 求圆和直线的两个交点。先建立方程组函数文件fxyz.m。function F=fxyz(X)x=X(1);y=X(2);z=X(3);F(1)=x2+y2+z2-9;F(2)=3*x+5*y+6*z;F(3)=x-3*y-6*z-1;再在MATLAB
7、命令行窗口，输入如下命令： X1=fsolve(fxyz,-1,1,-1,optimset(Display, off) %求第一个交点 X2=fsolve(fxyz,1,-1,1,optimset(Display, off) %求第二个交点例8-13 求函数在区间(-10，-1)和(1，10)上的最小值点。命令如下： f=(x) x-1./x+5; x,fmin=fminbnd(f,-10,-1) %求函数在(-10，-1)内的最小值点和最小值 fminbnd(f,1,10) %求函数在(1，10)内的最小值点例8-14 求函数f在(0.5，0.5，0.5)附近的最小值。建立函数文件fxyz0
8、.m。function f=fxyz0(u)x=u(1);y=u(2);z=u(3);f=x+y.2./x/4+z.2./y+2./z;在MALAB命令行窗口，输入如下命令： U,fmin=fminsearch(fxyz0,0.5,0.5,0.5) %求函数的最小值点和最小值例8-15 求解有约束最优化问题。首先编写定义目标函数函数文件fop.m。function f=fop(x)f=0.4*x(2)+x(1)2+x(2)2-x(1)*x(2)+1/30*x(1)3;再设定约束条件，并调用fmincon函数求解此约束最优化问题，程序如下：x0=0.5;0.5;A=-1,-0.5;-0.5,-1
9、;b=-0.4;-0.5;lb=0;0;option=optimset; option.LargeScale=off; option.Display =off;x,f=fmincon(fop,x0,A,b,lb,option)例8-16 求解线性规划问题。 命令如下：f=-5;-4;-6;A=1,-1,1;3,2,4;3,2,0;b=20;42;30;Aeq=;Beq=;LB=zeros(3,1);x,favl=linprog(f,A,b,Aeq,Beq,LB)例8-17 设有初值问题，试求其数值解，并与精确解相比较（精确解为y(t)=）。 建立函数文件funt.m。function yp=f
10、unt(t,y)yp=(y2-t-2)/4/(t+1);求解微分方程，程序如下：t0=0;tf=10;y0=2;t,y=ode23(funt,t0,tf,y0); %求数值解y1=sqrt(t+1)+1; %求精确解plot(t,y,b.,t,y1,r-) %通过图形来比较例8-18 已知一个二阶线性系统的微分方程为：取a=2，绘制系统的时间响应曲线和相平面图。函数ode23和ode45是对一阶常微分方程组设计的，因此对高阶常微分方程，需先将它转化为一阶常微分方程组，即状态方程。令x2=x，x1=x，则得到系统的状态方程：建立函数文件sys.m。function xdot=sys(t,x)xdot=-2*x(2);x(1); 取t0=0，tf=20，求微分方程的解，程序如下：t0=0;tf=20;t,x=ode45(sys,t0,tf,1,0); 程序运行后，查看结果。 t,x为直观地表示方程的解，可以绘制方程的时间响应曲线及相平面曲线（如图8-3所示），程序如下：subplot(2,2,1);plot(t,x(:,2) %系统时间响应曲线，即t-xsubplot(2,2,2);plot(x(:,
