计算方法大作业—质心型拉格朗日插值算法与Vandermonde矩阵逆矩阵算法

东南大学仪器科学与工程学院戴天宇1计算方法算法设计实践戴天宇22011229东南大学仪器科学与工程学院戴天宇2一、拉格朗日插值算法的改进1、质心型拉格朗日插值算法原理2、质心型拉格朗日插值算法讨论3、质心型拉格朗日插值算法Matlab实现4、数值算例二、Vandermonde矩阵逆矩阵的快速求法1-1、算法一原理阐述1-2、算法一Matlab实现1-3、算法一数值算例以及优缺点2-1、算法二原理阐述2-2、算法二Matlab实现2-3、算法二数值算例以及优缺点3、结论东南大学仪器科学与工程学院戴天宇3一、拉格朗日插值算法改进：1、原理阐释：𝜑𝑛(𝑥)=∏(𝑥−𝑥𝑗)𝑛𝑗=1东南大学仪器科学与工程学院戴天宇42、问题：1）、上述2种方法的工作量如何多少（如何尽可能的降低工作量）？解答：方式一：首先，计算Φn（x）需要n次乘法其次，每次计算wj需要n次乘法与加法再其次，乘n次，而后连加n次故，所需工作量为：（连乘n次+n次加法）*连加n次+乘法n次+连乘n次=2n^2+2n次即，2n^2量级的工作量方式二：首先，每次计算wj需要n次乘法与n次加法其次，将wj保存后，算分子与分母的加法各需n次，分子还需n次乘法故，所需工作量为：（连乘n次+n次加法）*连加n次+乘法n次+连加2n次=2n^2+3n次即，2n^2量级的工作量2）、上述两种方法的优缺点是什么？解答：方式一：优点是不需在添加节点时重新计算每个插值多项式，降低了一个量级的计算量，方式二：在方式一优点的基础上，只需要在计算wj的时候用到连乘，省去了计算Φn（x）的步骤，实现了wj的复用。并且实际上乘法比加法要更为复杂，所以只需要进行一次乘法的此方式更加优越。东南大学仪器科学与工程学院戴天宇53）、如果是均匀节点，wj应当如何算？解答：由于是均匀节点，故步长均一致，此时的计算公式大幅简化。令Δ为相邻节点差值，即为步长，则：Wj=1∏(𝑗−𝑖)𝛥𝑖≠𝑗=1(𝑗−1)𝛥∗(𝑗−2)𝛥∗…∗𝛥∗(−𝛥)∗(−2𝛥)…∗(j−n)𝛥=1(𝑗−1)∗(𝑗−2)∗…∗1∗(−1)∗(−2)…∗(j−n)∗𝛥𝑛−1首先令j=1，可以得到𝑊1=1(−1)∗(−2)…∗(1−n)∗𝛥𝑛−1若要计算W2，只需乘上一个因子，令其为V2：V2=1−𝑛2−1则：W2=W1∗V2故定义Vj为递推因子:Vj=(𝑗−1)−𝑛𝑗−1则：Wj=W1∗V2∗V3∗…∗Vj东南大学仪器科学与工程学院戴天宇6即，对于均匀节点，实际上计算所有的Wj只需要大概3n级别的运算量。又因为,在等距时,Wj是对称相反的的，即：𝑊1+𝑘=−𝑊𝑛−𝑘(k∈[1,n/2])故，实际上，只需要3*n/2级别的运算量3、编程序实现上述两种方法之一。实现算法二的程序。开始使用CPP进行编程，但结果由于CPP本身的精度问题，若要优化较为复杂，故后期转到MATLAB编程,但仍然具有输入点上限，质心型大概为五百个点，原算法更少，否则会溢出。于是，数值试验中输入点不会超过200个，以保证检验的正确。分别实现了原拉格朗日插值与算法二（质心拉格朗日插值），并作出了对比，MATLAB源代码如下：①、质心型拉格朗日插值：functionw=center_lagrange_w(x)%计算质心型拉格朗日插值算法中的w矩阵（wij）inum=length(x);w_n=ones(1,inum);%w矩阵的各元素取倒数后的倒数矩阵w=ones(1,inum);forcon_w=1:inumforcon_wj=1:(con_w-1)w_n(con_w)=w_n(con_w)*(x(con_w)-x(con_wj));endforcon_wj=(con_w+1):inumw_n(con_w)=w_n(con_w)*(x(con_w)-x(con_wj));endend%计算w_nforcon_wj=1:inumw(con_wj)=1/w_n(con_wj);%计算wendfunctionf_out=center_lagrange_out(x,f,w,x_in)%计算所求点函数值东南大学仪器科学与工程学院戴天宇7inum=length(x);onum=length(x_in);f_out=ones(1,onum);%输出矩阵初始化forcon=1:onumx_div=x_in(con)-x;div=w./x_div;sum1=sum(div.*f);sum2=sum(div);f_out(con)=sum1/sum2;%到此，f_out计算完毕forcon_x=1:inumifx_div(con_x)==0f_out(con)=f(con_x);%若输入值为插值点，输出原函数值endendendfunctionf_out=center_lagrange(x,f,x_in)%质心型拉格朗日插值w=center_lagrange_w(x);f_out=center_lagrange_out(x,f,w,x_in);end②、原拉格朗日插值算法：functionf=lagrange_old(x,y,x_in)%原拉格朗日插值函数onum=length(x_in);inum=length(x);f=zeros(1,onum);f=f.*0;%输出矩阵初始化fork=1:onumfori=1:inuml=y(i);forj=1:i-1l=l*(x_in(k)-x(j))/(x(i)-x(j));end;forj=i+1:inuml=l*(x_in(k)-x(j))/(x(i)-x(j));end;f(k)=f(k)+l;%插值、输出算法东南大学仪器科学与工程学院戴天宇8endend4、数值实验：实验一：输入数据：X=0:0.1：10（101个点）；f=sin(x);x_in=0:0.05:10(201个点);计算结果：时间：质心型——Elapsedtimeis0.018069seconds.原始型——Elapsedtimeis0.063726seconds.可见，点少时二者已有六倍速度差距。图像：质心型：原始型：东南大学仪器科学与工程学院戴天宇9可见，二者皆出现了龙格现象，这是由于插值点取值不够好所致，但是，质心型的两端相比原始型较为稳定。实验二：输入数据：X=0:0.05：10（201个点）；f=sin(x);x_in=0:0.0001:10(100001个点);计算结果：时间：质心型——Elapsedtimeis0.698409seconds.原始型——Elapsedtimeis112.357528seconds.此时，二者的速度差距已然差了许多个量级。证明了质心型插值算法的优越性。图像：质心型：原始型：东南大学仪器科学与工程学院戴天宇10与实验一一致，质心型数值稳定性更好。实验三：输入数据：X=0.05:0.05：10（200个点）；f=x.^(sin(x)+cos(x))/x+x.^2-log(x);x_in=0:0.001:10(10001个点);计算结果：时间：质心型——Elapsedtimeis0.080405seconds.原始型——Elapsedtimeis11.047469seconds.此时，二者的速度差距已然差了许多个量级。证明了质心型插值算法的优越性。图像（蓝色为原始图像，红色为插值多项式图像）：质心型原始型：东南大学仪器科学与工程学院戴天宇11实验四：输入数据：X=0:0.05：10（201个点）；f=exp(cosh(x)+sinc(x)-sinh(x)+sin(x));x_in=0:0.001:10(10001个点);计算结果：时间：质心型——Elapsedtimeis0.160583seconds.原始型——Elapsedtimeis11.047469seconds.此时，二者的速度差距已然差了许多个量级。证明了质心型插值算法的优越性。图像（蓝色为原始图像，红色为插值多项式图像）：质心型：原始型：以上实验均说明，质心型的速度与数值稳定性都比原始型好非常之多。东南大学仪器科学与工程学院戴天宇12二、范德蒙德矩阵的逆矩阵的快速算法：构造一个算法，使得即，在6n^2的乘除法运算内算出一个标准范德蒙德矩阵的逆矩阵。解答：东南大学仪器科学与工程学院戴天宇131、算法一：1)、原理阐释：参考陆全、任学明的《一类广义范德蒙矩阵求逆的快速算法》，得出以下解决方法：有：此例中，标准范德蒙德矩阵的定义为：𝑉𝑠0(𝑛)=[1𝑐1c121𝑐2c221𝑐3c32⋯c1𝑛−3c1𝑛−2c1𝑛c2𝑛−3c2𝑛−2c2𝑛c3𝑛−3c3𝑛−2c3𝑛⋮⋱⋮1𝑐𝑛−2c𝑛−221𝑐𝑛−1c𝑛−121𝑐𝑛c𝑛2⋯c𝑛−2𝑛−3c𝑛−2𝑛−2c𝑛−2𝑛c𝑛−1𝑛−3c𝑛−1𝑛−2c𝑛−1𝑛c𝑛𝑛−3c𝑛𝑛−2c𝑛𝑛]𝑇记v𝑗为𝑉𝑠0(𝑛)−1的第0个列向量。则有递推公式：{𝑣1=𝑤𝑣𝑗=𝐷∗𝑣𝑛=𝑢𝑣𝑗−1−𝑥𝑗−1∗𝑤−𝑦𝑗−1∗𝑢(𝑗=2,…𝑛−1)其中：D=diag(1𝑎1,1𝑎2,…,1𝑎𝑛)而计算w、u、x、y的算法如下：东南大学仪器科学与工程学院戴天宇142)、代码实现：function[V_inv,amount]=Vandermonde_inv_fast(V)n=length(V);%获取待求矩阵长度C=V(2,:);%c1,c2,...,cnD=diag(1./C);%构造对角矩阵DV_inv=ones(n,n);%逆矩阵初始化amount=0;%乘除法次数计数初始化w=zeros(n,n);u=zeros(n,n);x=zeros(n,n);y=zeros(n,n);w(1,1)=1;x(1,1)=1/C(1);y(1,1)=C(1)^(n-1)-C(1)^(n-2);%wuxy初始化东南大学仪器科学与工程学院戴天宇15l=zeros(1,n);d=zeros(1,n);t=zeros(1,n);p=zeros(n,n);q=zeros(n,n);%ldtpq初始化%wuxy计算开始fork=2:(n-1)forj=1:(k-1)l(k)=l(k)+C(j)^(k-1)*w(k-1,j);d(k)=d(k)+C(k)^(j-1)*x(k-1,j);t(k)=t(k)+C(k)^(j-1)*y(k-1,j);amount=amount+3;endl(k)=-l(k);d(k)=1/C(k)-d(k);t(k)=C(k)^(n-1)-C(k)^(n-2)-t(k);p(k,k)=1/(C(k)*l(k)*d(k));amount=amount+4;q(k,k)=p(k,k);fori=1:(k-1)p(k,i)=p(k,k)*d(k)*w(k-1,i)/(1/C(i)-1/C(k));amount=amount+2;endfori=2:(k-1)q(k,i)=-q(k,k)*l(k)*x(k-1,i-1);amount=amount+2;endforj=2:kq(k,1)=q(k,1)+C(1)^(j-1)*q(k,j);amount=amount+1;endq(k,1)=-q(k,1);w(k,:)=w(k-1,:)+l(k)*p(k,:);x(k,:)=x(k-1,:)+d(k)*q(k,:);y(k,:)=y(k-1,:)+t(k)*q(k,:);amount=amount+3*k;endforj=1:(n-1)l(n)=l(n)+C(j)^n*w(n-1,j);d(n)=d(n)+C(n)^(j-1)*x(n-1,j);东南大学仪器科学与工程学院戴天宇16t(n)=t(n)+C(n)^(j-1)*y(n-1,j);amount=amount+3;endl(n)=-l(n);d(n)=1/C(n)-d(n);t(n)=C(n)^(n-1)-C(