材料0901班杨东麟0604090113(Page41)

１、编程实现以下科学计算算法，并举一例应用之。

（参考书籍《精 通ＭＡＬＡＢ科学计算》 ，王正林等著，电子工业出版社，２００９ 年）“Gauss-Seidel 迭代法线性方程组求解”（1）迭代解法的基本思想：根据给定方程组，设计出一个迭代公式，构造一数组的序列xi 0,代入迭代公式，计算出xi 1,在代入迭代公式，经过k 次迭代运算后得到xi k ,若xi k 收敛于某一极限数组xi ，则xi 就是方程组的近似解。

迭代过程本质上就是计算极限的过程，一般不能得到精确解。

但迭代的优点是程序简单，适合大型方程组求解，然而，缺点是要判断迭代是否收敛和收敛速度的问题。

(2)算法说明：Gauss-Seidel 迭代法与简单迭代法类似，只是迭代公式有所改进。

简单迭代法：kjnij j iiij iii k ix A A A B x∑≠=+-=,11；Gauss-Seidel 迭代法：kjni j iiij k ji ij iiij iii k ix A A xA A AB x∑∑+=+-=+--=1111；设方程组Ax ＝b ，其中A 和b 中的元素都为常数，且A 为非奇异，则A 分可写成：A=D-L-U 。

其中D 上网意义同Jacobi 迭代法，L 为下三角矩阵，U 为上三角矩阵，他的迭代公式为：bL D Ux L D x k k 111)()(--+-+-=在MA TLAB 中编程实现的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%求出所需精度的解实际的迭代步数：nif nargin ==3eps=1.0e-6; %eps表示迭代精度M=200; %M表示迭代步数的限制值elseif nargin ==4M=200;elseif nargin <3errorreturn;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;%迭代过程while norm(x-x0)>epsx0=x;x=G*x0+f;n=n+1; %n为最终求出解时的迭代步数if n>=Mdisp('Warning:迭代次数太多，可能不收敛！');endend(4)进行实例分析：>>A=[1.4449 0.7948 0.8801;0.6946 1.9568 0.1730;0.6213 0.5226 1.9797]; >> b=[1 0 1]';>> 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） 当t<0时，开关K 位于“1”，电路以达到平衡U c （0+）=Uc （0-）=-12V ， i R2（0-）=3A ； （2） 当t>0时，i c （0+）=-U c （0+）/(R 2*R 3/（R 2+R 3）)，达到稳态后，电容中将无电流i cf =0A ，电流源的全部电流将在两个电阻之间分配，保证端电压相同(也就是电容上的终电压)。

即：U R3= U cf =12V 。

（3） 由一阶响应电路可用三要素法得到电压公式：τtf f eu u u t u -+-+=))0(()(，时间常数RC =τ。

τtcf c cf c eU U U t U -+-+=))0(()(，且R eq =12*6/18=4，RC =τ=4(s)。

同时，以相同的方法得到电流公式：τtf c c f c c e i i i t i -+-+=))0(()(，RC =τ=4(s)。

（4） 可以得到电阻R 2的表达式：[]3232)()(R R R t i I t I c s R +⋅-=。

二、源程序设计：>> 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...)(a x a xa x a x y n n n n ++++=--设有n+1个互异节点n x x x <<<...10，且)(i i x f y =，i 为整数，构造L n (x),使L n (x j )=y j ，j 为整数。

定义，若n 次多项式l j (x)在n+1个节点n x x x <<<...10上满足条件⎩⎨⎧=01)(j k x l ，为1时，j=k;为0时，j~=k;则称这n+1个n 次多项式l 0(x),l 1(x),…,l n (x)为节点x 0,x 1,…,x n 上的n 次插值基函数。

即为：∏≠=--=nki i i ki k x xx x x l 0)()()(，k=(0,1,2,…,n);故满足插值条件多项式为：n n n y x l y x l y x l x L )(...)()()(1100+++=，称其为Lagrange 插值多项式。

○2n-1阶拉格朗日插值matlab 源程序： function y=lagrange(x0,y0,x) % 计算n 次lagrange 多项式插值 % x0,y0输入插值节点向量 % x 以向量形式输入的插值点 % y 输出的插值点函数数值 n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=kp=p*(z-x0(j))/(x0(k)-x0(j)); end ends=p*y0(k)+s; end y(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?输出：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流程图：开始输入：等距插值点步长d=2*pi/8输入：x,y=sin(x)得到插值点用Lagrange进行八阶九点差值计算绘图比较sin(x)与Lagrange插值图像比较两图误差结束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 =0 1.2500 2.5000 3.7500 5.0000y =1.0000 0.2748 0.1414 0.0937 0.0698（3）Lagrange差值多项式：L n(x)=l0(x)*y0+ l1(x)*y1+ l2(x)*y2+ l3(x)*y3+ l4(x)*y4（4）间隔0.2计算的误差值：e =Columns 1 through 90 -0.0209 -0.0405 -0.0406 -0.0294 -0.0152 -0.0026 0.0063 0.0111Columns 10 through 180.0122 0.0104 0.0067 0.0023 -0.0021 -0.0057 -0.0078 -0.0082 -0.0066Columns 19 through 26-0.0034 0.0012 0.0065 0.0115 0.0150 0.0155 0.0112 0 （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的误差分布图结束。