MATLAB及其在电气工程中的应用 作者 苏小林 第三章.pptVIP

第3章 数值计算与分析 课堂练习： 1、有一电路的节点电压方程如下，求该电路的各节点电压。 2、下面是某物理量y在几个不同时间点的观察值记录： t y 0.0 0.82 0.3 0.72 0.8 0.63 1.1 0.60 1.6 0.55 2.3 0.50 用衰减指数函数y=c1+c2e-t建立这组数据的模型。 3、对矩阵进行奇异值分解。 A=[1 3;2 4;5 7;6 8]; 课堂练习：1、求矩阵的特征值和特征向量。 A=[2 -4 -1;4 1 -10;-3 15 -10] 在初始条件 ， ，且 情况下的解，并用图示显示解的结果。 (1) 把高阶微分方程改写成一阶微分方程组 令 ， ，原二阶微分方程改写为如下一阶微分方程组： ， (2) 编写一阶微分方程组的函数文件 function dydt=vdp1(t,y) mu=1; dydt=[y(2);mu*(1-y(1)^2)*y(2)-y(1)]; 将上面的程序代码保存为vdp1.m。 (3) 根据微分方程，确定需采用的算法，对于范德波耳非刚性方程，选择ode45解算微分方程，解算时间间隔取[0,20] [t,y]=ode45(@vdp1,[0,20],[2;0]); 数组y的第1列为y1，第2列为y2，每行为它们在对应时间点的解。解法函数的第1个输入参数为定义微分方程的函数句柄。 (4) 用图示观察解的结果 plot(t,y(:,1),-,t,y(:,2), --) title( Solution of van der Pol Equation,\mu=1 ) xlabel( time t ) ylabel( solution y ) legend( y_1 , y_2 ) 图3-2 微分方程的解 【例3-30】将例3-29的范德波耳微分方程中的m取为1000时，微分方程成为刚性的。求该刚性范德波耳微分方程的解。 (1) 编写一阶微分方程组的函数文件vdp1000 function dydt=vdp1000(t,y) mu=1000; dydt=[y(2);mu*(1-y(1)^2)*y(2)-y(1)]; (2) 针对刚性范德波耳微分方程，采用ode15s解法，时间间隔为 [0,3000] tspan=[0 3000]; y0=[2;0]; [t,y]=ode15s(@vdp1,tspan,y0); 课堂练习： 1、求解积分 的数值。 2、求解二重积分 的数值。 3、求解范德波耳(van der Pol)方程 在初始条件 且 情况下的解，并用图示显示解的结果。 4、将上题的范德波耳微分方程中的取为1000时，微分方程成为刚性的。求该刚性范德波耳微分方程的解。 x=A\b % 求欠定方程组的解 x = -0.2798 0.4940 0 0.4345 3.2.2 矩阵分解 在求解线性方程组、进行模型参数估计等许多问题中，常需要将矩阵分解。所谓矩阵分解主要是指将一个矩阵分解为几个比较简单的矩阵连乘的形式。在MATLAB中，线性方程组的求解主要基于三种基本的矩阵分解：Cholesky分解、LU分解和QR分解。对于这些分解，MATLAB都提供了对应的函数来实现。此外，奇异值分解在在矩阵分析中也有着非常重要的作用。 1. Cholesky分解 对于正定对称矩阵A，可以分解为一个上三角矩阵R和其转置矩阵的乘积，即A=RTR，这种分解称为Cholesky分解。在MATLAB中，使用函数chol来实现。其常用命令格式为R=chol(A)，其中A是对称的正定矩阵，R是求得的上三角矩阵；若矩阵A是非正定矩阵，该命令返回错误信息。 对于线性方程组Ax=b，若A是对称的正定矩阵，其A=RTR，则线性方程组可以转换为RTRx=b，该方程组的解为x=R\(RT\b)，MATLAB命令为x=R\(R’\b)。它与左除x=A\b所得结果完全相同，但对于n×n的矩阵A，命令ch