材料090班杨东麟0604090113(Page41)

１、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“Gauss-Seidel迭代法线性方程组求解”（1）迭代解法的基本思想：根据给定方程组，设计出一个迭代公式，构造一数组的序列xi0,代入迭代公式，计算出xi1,在代入迭代公式，经过k次迭代运算后得到xik,若xik收敛于某一极限数组xi，则xi就是方程组的近似解。迭代过程本质上就是计算极限的过程，一般不能得到精确解。但迭代的优点是程序简单，适合大型方程组求解，然而，缺点是要判断迭代是否收敛和收敛速度的问题。(2)算法说明：Gauss-Seidel迭代法与简单迭代法类似，只是迭代公式有所改进。简单迭代法：kjnijjiiijiiikixAAABx,11；Gauss-Seidel迭代法：kjnijiiijkjiijiiijiiikixAAxAAABx1111；设方程组Ax＝b，其中A和b中的元素都为常数，且A为非奇异，则A分可写成：A=D-L-U。其中D上网意义同Jacobi迭代法，L为下三角矩阵，U为上三角矩阵，他的迭代公式为：bLDUxLDxkk111)()(在MATLAB中编程实现的Gauss-Seidel迭代法函数为：gauseidel。功能：用Gauss-Seidel迭代求线性方程组ax=b的解。调用格式：[x,n]=gauseidel(A,b,x0,eps,M).其中，A为线性方程组的系数矩阵；b为线性方程组中的常数向量；x0为迭代初始向量；eps为解的精度控制（此参数可选）；M为迭代步数控制（此参数可选）；x为线性方程组的解；n为求出所需精度的解实际迭代步数。(3)Gauss-Seidel迭代法的MATLAB程序代码如下:function[x,n]=gauseidel(A,b,x0,eps,M)%采用Gauss-Seidel迭代法求线性方程组Ax=b的解%线性方程组的系数矩阵：A%线性方程组的常数向量：b%迭代初始向量：x0%解的精度控制：eps%迭代步数控制：M%线性方程组的解：x%求出所需精度的解实际的迭代步数：nifnargin==3eps=1.0e-6;%eps表示迭代精度M=200;%M表示迭代步数的限制值elseifnargin==4M=200;elseifnargin3errorreturn;endD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;%迭代过程whilenorm(x-x0)epsx0=x;x=G*x0+f;n=n+1;%n为最终求出解时的迭代步数ifn=Mdisp('Warning:迭代次数太多，可能不收敛！');endend(4)进行实例分析：A=[1.44490.79480.8801;0.69461.95680.1730;0.62130.52261.9797];b=[101]';x0=zeros(3,1);[x,n]=gauseidel(A,b,x0)x=%输出结果0.5929-0.24440.3836n=11%输出迭代次数n(5)运行图以即流程图：①运行图：②流程图：x=G*x0+fnorm(x-x0)eps？n=200？Warning:迭代次数太多，可能不收敛x0=x;n=n+1最大步数M为200是开始结束读取数据nargin==3？输出结果n=1eps=10--6，最大步数M=200否否是是否②例题流程图：输入系数矩阵A输入常数向量b和初始向量x0[x,n]=gauseidel(A,b,x0)输出计算结果输出迭代次数一、分析电路：（RC充电电路）（1）当t0时，开关K位于“1”，电路以达到平衡Uc（0+）=Uc（0-）=-12V，iR2（0-）=3A；（2）当t0时，ic（0+）=-Uc（0+）/(R2*R3/（R2+R3）)，达到稳态后，电容中将无电流icf=0A，电流源的全部电流将在两个电阻之间分配，保证端电压相同(也就是电容上的终电压)。即：UR3=Ucf=12V。（3）由一阶响应电路可用三要素法得到电压公式：tffeuuutu))0(()(，时间常数RC。tcfccfceUUUtU))0(()(，且Req=12*6/18=4，RC=4(s)。同时，以相同的方法得到电流公式：tfccfcceiiiti))0(()(，RC=4(s)。（4）可以得到电阻R2的表达式：3232)()(RRRtiItIcsR。二、源程序设计：r1=3;r2=12;r3=6;us=18;is=3;c=1;uc0=-12;%电容C的初始电压值ucf=12;%电容C的最终稳态电压值icf=0;%电容C的最终稳态电电流值T=r2*r3/(r2+r3);%电路时间常数t=[0:0.01:25];%时间变量uct=ucf+(uc0-ucf)*exp(-t./T);%计算电容C的响应电压ic0=-uc0/(r2*r3/(r2+r3));%电容C的初始电流值ict=icf+(ic0-icf)*exp(-t./T);%计算电容C的响应电流ir20=3;%电阻R2的初始稳态电流ir2t=(is-ict)*r3/(r2+r3);%计算电阻R2的响应电流subplot(3,1,1);plot(t,uct,'-r');title('Uc(t)');subplot(3,1,2);plot(t,ict,'-g');title('Ic(t)');subplot(3,1,3);plot(t,ir2t,'-b');title('Ir2(t)');%输出Uc（t）和IR2（t）的响应图像三、波形图：四、流程图：开始读取已知数据计算稳态电流电压数据运用响应电路三要素构造函数uct=ucf+(uc0-ucf)*exp(-t./T)、ict=icf+(ic0-icf)*exp(-t./T)和ir2t=(is-ict)*r3/(r2+r3)结束输出Uc（t）、Ic(t)和IR2（t）的响应图像1.分别用五点和九点等距插值：○1Lagrange插值多项式原理：n次插值就是利用n+1个插值节点构造的n次的插值多项式，根据差值函数条件，可以得到：01111...)(axaxaxaxynnnn设有n+1个互异节点nxxx...10，且)(iixfy，i为整数，构造Ln(x),使Ln(xj)=yj，j为整数。定义，若n次多项式lj(x)在n+1个节点nxxx...10上满足条件01)(jkxl，为1时，j=k;为0时，j~=k;则称这n+1个n次多项式l0(x),l1(x),…,ln(x)为节点x0,x1,…,xn上的n次插值基函数。即为：nkiiikikxxxxxl0)()()(，k=(0,1,2,…,n);故满足插值条件多项式为：nnnyxlyxlyxlxL)(...)()()(1100，称其为Lagrange插值多项式。○2n-1阶拉格朗日插值matlab源程序：functiony=lagrange(x0,y0,x)%计算n次lagrange多项式插值%x0,y0输入插值节点向量%x以向量形式输入的插值点%y输出的插值点函数数值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○3流程图：（1）五点等距插值的四阶Lagrange多项式○1源程序：d=2*pi/4;x=[0:d:2*pi];y=sin(x);x0=[0:0.1:2*pi];y0=lagrange(x,y,x0);y1=sin(x0);输入：z输入：x0,y0开始s=0k=1p=1j=1p=p*(z-x0(j))/(x0(k)-x0(j))(j~=k)j=j+1j=n?s=p*y0(k)+sk=n?=n?输出：z,s结束k=k+1是是否否plot(x0,y1,'--r',x0,y0,'-b',x,y,'o')legend('sin(x)','四阶插值多项式','lagrange插值点')○2运行图：○3流程图：输入：x,y=sin(x)得到插值点用Lagrange进行四阶五点差值计算输入：等距插值点步长d=2*pi/4开始绘图比较sin(x)与Lagrange插值图像比较两图误差结束（2）九点等距插值的四阶Lagrange多项式○1源程序：d=2*pi/8;x=[0:d:2*pi];y=sin(x);x0=[0:0.1:2*pi];y0=lagrange(x,y,x0);y1=sin(x0);plot(x0,y1,'-r',x0,y0,'-b',x,y,'+')legend('sin(x)','八阶插值多项式','lagrange插值点')○2运行图：○3流程图：2.四阶Lagrange插值多项式近似：（1）源程序：d=5/4;x=[0:d:5];y=(1+x)./(1+2*x+3*x.^2);x0=[0:0.05:5];y0=lagrange(x,y,x0);y1=(1+x0)./(1+2*x0+3*x0.^2);x1=[0:0.2:5];y2=(1+x1)./(1+2*x1+3*x1.^2);y3=lagrange(x,y,x1);e=y2-y3;plot(x,y,'o',x0,y1,'-r',x0,y0,'-b',x1,e,'+')legend('插值点','原函数图象','四阶Lagrange插值多项式','误差分布图')（2）插值点：x=01.25002.50003.75005.0000结束绘图比较sin(x)与Lagrange插值图像比较两图误差输入：等距插值点步长d=2*pi/8输入：x,y=sin(x)得到插值点用Lagrange进行八阶九点差值计算开始y=1.00000.27480.14140.09370.0698（3）Lagrange差值多项式：Ln(x)=l0(x)*y0+l1(x)*y1+l2(x)*y2+l3(x)*y3+l4(x)*y4（4）间隔0.2计算的误差值：e=Columns1through90-0.0209-0.0405-0.0406-0.0294-0.0152-0.00260.00630.0111Columns10through180.01220.01040.00670.0023-0.0021-0.0057-0.0078-0.0082-0.0066Columns19through26-0.00340.00120.00650.01150.01500.01550.01120（5）函数图象以及误差分布图：（6）流程图：开始输入：等距插值点步长d=5/4输入：x,y=(1+x)./(1+2*x+3*x.^2)确定插值点确定Lagrange多项式的y(i)值,i=0~4，即得到多项式用Lagrange进行四阶5点差值计算绘图比较Lagrange插值多项式图像与原函数图y2=(1+x1)./(1+2*x1+3*x1.^2)以0.2为步长，计算每隔0.2的点的误差值e=y2-y3绘出e=y2-y3的误差分布图结束