矩阵与数值分析

矩阵与数值分析

学 院 专 业 班 级 学 号 姓 名

电子信息与电气工程学部

生物医学工程

刘江涛

1：考虑计算给定向量的范数；输入向量x?(x1,x2,?,xn)T，输出x1,x2,x?，请编制一个通用程序，并用你编制的程序计算如下向量的范数：

1??11Tx??1,,,?,?,y??1,2,?,n?

n??23对n?10,100,1000甚至更大的n计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？

通用求范数程序： function NORM(x) y1=sum(abs(x)); y2=(sum(x.^2))^(1/2); y3=max(abs(x));

fprintf('1-范数=%g； 2-范数= %g； inf-范数=%g\\n',y1,y2,y3); 例题的运行程序： function xianglaing(n) x=[]; y=[]; for i=1:n x(i)=1/i; y(i)=i; end

disp('x 的范数：'); NORM(x'); disp(' ')

disp('y 的范数：'); NORM(y'); 运行结果如下表：

T n 范数 10 100 1000 x1 2.92897 5.18738 7.48547 x2 x? 1 1 1 1.2449 1.27866 1.28216 n 范数 10 100 1000 y1 55 5050 500500 y2 y? 10 100 1000 19.6214 581.679 18271.1 根据上述的两个表的运行结果，我们可以得知无论n的值如何变化，对于x??1恒成立；y??n恒成立，其1-范数与2-范数随着n的增大而增大，但是其变化越来越小，这是因为计算在进行数值计算时有误差存在，对于表达式(1)当n很大时

1却很n小，会出现“大数吃小数的现象”；修改方案：当n很大时我们避免用n做除数，因为当n非常大时

1?0成立；所以在求解其范数时我们从小数开始相加，无穷个非常n小的数值相加也可能是个很大的数，从而可以避免两个数相加时出现“大数吃小数”的现象；

2：考虑y?f(x)?ln(1?x)，其中定义f(0)?1，此时f(x)是连续函数，用此公x式计算当x?[?10?15,10?15]时的函数值，画出图像。另一方面，考虑下面算法：

d?1?x;ifd?1theny?1elsey?lnd/(d?1)endif

用此算法计算x?[?10?15,10?15]时的函数值，画出图像，比较一下发生了什么？ 程序：

x=-10^(-15):10^(-20):10^(-15);

if (x==0) f=1; else f=log(1+x)/x; end figure(1) plot(x,f); d=1+x; if d==1 y=1; else

y=log(d)/(d-1); end figure(2) Plot(x,y);

有图可知，直接用公式f(x)?ln(1?x)计算x?[?10?15,10?15]的函数值时，除了在xx?0出的值为1，其他的值都是无限趋近于1；而利用算法二算出的结果全为1；出

现这这情况的原因是x的取值非常接近于0，在用公式d?1?x求d得过程中出现了大数吃小数的情况，所以在用计算机计算时d?1恒成立，从而使y?1恒成立；

3： 首先编写一个利用秦九韶算法计算一个多项式在定点的函数值的通用程序，你的程序包括输入多项式的系数以及定点，输出函数值，利用你编写的程序计算

f(x)?(x?2)9?x9?18x8?144x7?672x6?2016x5?4032x4?5376x3?4608x2?2304x?512在x?2邻域附近的值，画出p(x)在x?[1.95,20.5]上地图像。

秦九韶算法的通用程序：

%A为多项式的以升幂排列的系数，x为初始值 function p=qinjiushao(A,x) a=A; [~,n]=size(a); n=n-1; S=[];

S(n+1)=a(n+1); for k=n:-1:1

S(k)=x.*S(k+1)+a(k);

end p=S(1);

利用上述程序计算p(x)在x?2邻域附近的值具体见下表：

x p(x) 1.95 -9.6634e-13 1.96 -1.1255e-11 1.97 2.8422e-12 1.98 7.3896e-12 1.99 5.3433e-12 2.00 0 x p(x) 2.01 -3.7517e-12 2.02 -7.7875e-12 2.03 4.4338e-12 2.04 -3.6380e-12 2.05 5.1159e-12 当x?[1.95,20.5]时，p(x)的图像如下： 画图程序如下： function huatu(A,x) [~,n]=size(x); for i=1:n

y(i)=qinjiushao(A,x(i)); end plot(x,y); 程序运行如下： >> x=1.95:0.01:20.5;

>> A=[-512 2304 -4608 5376 -4032 2016 -672 144 -18 1]; >> huatu(A,x)

3x 10112.521.5p(x)10.50-0.50510x152025

4：编制计算机给定矩阵A的LU分解和PLU分解的通用程序，然后利用你编写的程序完成下面两个计算任务：

考虑

?1??1?A??????1???10?1?????01??Rn?n

?11??11??0??????1??1自己取定x?Rn，并计算b?Ax。然后用你编制的不选主元的Gauss消去法求解

?。对n从5到30估计计算解的精度。 该方程组，记你计算出的解为x（2）对n从5到30计算出其逆矩阵。 LU分解的通用程序： function LU(A) [m,n]=size(A); L=zeros(m,n); U=zeros(m,n); for j=1:n

U(1,j)=A(1,j); L(j,j)=1;

end for j=2:m

L(j,1)=A(j,1)/U(1,1); end for i=2:n for j=i:n sum=0; for k=1:i-1

sum=sum+L(i,k)*U(k,j); end

U(i,j)=A(i,j)-sum; end for j=i+1:n sum=0; for k=1:i-1

sum=sum+L(j,k)*U(k,i); end

L(j,i)=(A(j,i)-sum)/U(i,i); end end disp('L ='); disp(L); disp('U ='); disp(U);

PLU分解的通用程序： function PLU(A)

[~,n]=size(A); Ip=1:n; for k=1:n-1

[~,r]=max(abs(A(k:n,k))); r=r+(k-1); if r>k

A([k,r],:)=A([r,k],:); Ip([k,r])=Ip([r,k]); end for p=k+1:n mu=A(p,k)/A(k,k); A(p,k)=mu;

A(p,k+1:n)=A(p,k+1:n)-mu*A(k,k+1:n); end end p=eye(n,n); P=zeros(n,n); for i=1:n

P(i,1:n)=p(Ip(i),1:n); end

L=tril(A,-1)+eye(n) ; U=triu(A); disp('P=') disp(P); disp('L=') disp(L);

disp('U='); disp(U);

我选取b?[123?n]Tn?[5,30],n?N*，则我们可以计算出其精确解

2n?1?12n?2?1?1x?[?n?1?n?2?2222n?1T]2n?1n?[5,30]n?N*

实现求解的程序如下： function INV(n) A=ones(n); for i=2:n for j=1:i-1

A(i,j)=-1*A(i,j); end for j=i:n-1 A(i-1,j)=0; end end for i=1:n b(i)=i;

X(i)=-(2^(n-i)-1)/(2^(n-i)); end

X(n)=(2^n-1)/(2^(n-1)); [L,U]=LU(A); x1=inv(U)*inv(L)*b' disp(x1);

我们利用上述的程序计算出的结果如下表：

n x1,x2,x3,?,xn

5 6 7 … -0.9375,-0.8750,-0.7500,-0.5000,1.9375 -0.9688,-0.9375,-0.8750,-0.7500,-0.5000,1.9688 -0.9844,-0.9688,-0.9375,-0.8750,-0.7500,-0.5000,1.9844 … 利用我们求出的精确解与用程序求出的近似解求其误差，再利用matlab编程实现时其误差为零。 （2）

对于题目中的A的求逆的通用程序： function INV(n) A=ones(n); for i=2:n for j=1:i-1

A(i,j)=-1*A(i,j); end for j=i:n-1 A(i-1,j)=0; end end B=inv(A);

disp('A 的逆为：'); disp(B); n=5时：

-0.2500 -0.1250 -0.0625 -0.0625?? 0.5000? ?0 0.5000 -0.2500 -0.1250 -0.1250???1? A?? 0 0 0.5000 -0.2500 -0.2500?? 0 0 0 0.5000 -0.5000???? 0.2500 0.1250 0.0625 0.0625? 0.5000?n=6时：

-0.2500 -0.1250 -0.0625 -0.0313 -0.0313?? 0.5000? 0 ?0.5000 -0.2500 -0.1250 -0.0625 -0.0625??? 0 ?0 0.5000 -0.2500 -0.1250 -0.1250A?1???

0 0 0.5000 -0.2500 -0.2500? 0 ?? 0 ?0 0 0 0.5000 -0.5000?? 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313????n=7时：

-0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156?? 0.5000? 0 ? 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313??? 0 ? 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625??A?1?? 0 0 0 0.5000 -0.2500 -0.1250 -0.1250?

? 0 0 0 0 0.5000 -0.2500 -0.2500??? 0 0 0 0 0 0.5000 -0.5000??? 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156???n=8时：

-0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0078 -0.0078?? 0.5000? 0 ? 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156??? 0 ? 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313?? 0 0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625? A?1??? 0 0 0 0 0.5000 -0.2500 -0.1250 -0.1250??? 0 0 0 0 0 0.5000 -0.2500 -0.2500??? 0 0 0 0 0 0 0.5000 -0.5000??? 0.2500 0.1250 0.0625 0.0313 0.0156 0.0078 0.0078??? 0.5000?.......

在此不再一一列举，都可以用上述程序算出；

5：编制计算对称正定阵的Cholesky分解的通用程序，并利用你编制的程序计算

Ax?b，其中A?(aij)?Rn?n,aij?1，b可以有你自己取定，对n从10到20验

i?j?1证程序的可靠性。

Cholesky分解求L的通用程序： %LT代表L的转置 function [L,LT]=Cholesky(A) [n,m]=size(A); L=zeros(n,n); for j=1:m sum=0; for k=1:j-1

sum=sum+(L(j,k))^2; end

L(j,j)=(A(j,j)-sum)^(1/2); for i=j+1:n sum=0; for k=1:j-1

sum=sum+L(i,k)*L(j,k); end

L(i,j)=(A(i,j)-sum)/L(j,j); end end L=L; LT=L';

求解Ax?b的方程解的程序如下： function cholesky_qiu_jie(n,b) A=zeros(n); for i=1:n for j=1:n

A(i,j)=1/(i+j-1); end end

[L,LT]=Cholesky(A); Y=inv(LT)*inv(L)*b; disp(Y'); >> n=10;

>> b=[1 0 0 0 0 0 0 0 0 0]'; >> cholesky_qiu_jie(n,b) 1.0e+06 *

Columns 1 through 5

0.000099995224850

-0.004949584178747

0.079191090719606

-0.600518630835813 2.522130446082295

Columns 6 through 10

-6.305225874989181

9.607833253959356

-8.749889022961289

4.374900125373453 -0.923581794866949

>> n=11;

>> b=[1 0 0 0 0 0 0 0 0 0 0]'; >> cholesky_qiu_jie(n,b) 1.0e+07 * Columns 1 through 5

0.000012085421839 -0.000724460241949 0.014116768083723

-0.131679518851255 0.690983916160226

Columns 6 through 10 -2.210252357004733

4.471585077358848

-5.747468746188864

4.548898816586863 -2.021271620300206

Column 11

0.385801129112326 n=12;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1]'; >> cholesky_qiu_jie(n,b) 1.0e+08 * Columns 1 through 5 -0.000000106586186

0.000015496760612

-0.000549242897797

0.008320106358643 -0.067089571289062

Columns 6 through 10 0.321428166562500

-0.969534469375000

1.888382749062500

-2.369824064062500 1.849525708125000

Columns 11 through 12

-0.816237543125000 0.155564209082031 n=13;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1]'; >> cholesky_qiu_jie(n,b) 1.0e+10 * Columns 1 through 5 0.000000017255662

-0.000002801898660

0.000111021739307

-0.001891219259375 0.017318444275000

Columns 6 through 10

-0.095564468900000 0.338540141400000 -0.795937993200000

1.255304630000000 -1.312819645600000

Columns 11 through 13

0.873192546400000 -0.334354244200000 0.056103628050000 >> n=14;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0]'; >> cholesky_qiu_jie(n,b) 1.0e+15 * Columns 1 through 5 -0.000000018973466

0.000003062445902

-0.000120619605275

0.002041840742992 -0.018571410648157

Columns 6 through 10 0.101712531234567

-0.357252637738148

0.831474891793232

-1.294880899877584 1.331275372322000

Columns 11 through 14 -0.862636505012200 -0.002171662407011

>> n=15;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0 1]'; >> cholesky_qiu_jie(n,b) 1.0e+16 * Columns 1 through 5 0.000000009159679

-0.000000881206821

0.000013799271738

0.314323195943382

-0.045197162809621

0.000103167625626 -0.003987948970647

Columns 6 through 10 0.038905667203020

-0.198056922710865

0.601726021228474

-1.119920061974954 1.191735365493764

Columns 11 through 15 -0.450794442163960

-0.496191645217859

0.756621451313167

-0.399838355711562 0.079684802839153

>> n=16;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0]'; >> cholesky_qiu_jie(n,b) 1.0e+16 * Columns 1 through 5 0.000000106506398

-0.000014879116455

0.000507146024044

-0.007257871409944 0.052560866946652

Columns 6 through 10 -0.198284962679135

0.290609837891524

0.581225032778005

-3.497134240175643 6.758520114521732

Columns 11 through 15 -5.326111865614911

-2.146668037387352

8.700448141028517

-8.161919041143481 3.588648336729691

Column 16

-0.635128701903556 >> n=17;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1]'; >> cholesky_qiu_jie(n,b) 1.0e+17 * Columns 1 through 5 0.000000014777084

-0.000002118394303

0.000074232548837

-0.001097606950101 0.008314003699029

Columns 6 through 10 -0.034029369522343

0.065917257262332

0.015959536689882

-0.377621602904773 0.819475115416494

Columns 11 through 15 -0.660163184213223

-0.295786673257216

1.100782605917994

-0.952907140993597 0.339855611203294

Columns 16 through 17

-0.014551513983896 -0.014219174116106 >> n=18;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 2]'; >> cholesky_qiu_jie(n,b) 1.0e+16 * Columns 1 through 5 0.000000155967115

-0.000021645380087

0.000727005315139

-0.010130391409408 0.069910333639844

Columns 6 through 10 -0.237692894985775

0.207119320267744

1.320710222270083

-5.362833277085945 8.761578864695567

Columns 11 through 15 -5.416385894431453

-3.678202876264916

9.162910938126201

-7.610028175492047 4.801505077240169

Columns 16 through 18

-3.460280780403269 1.887745959501761 -0.436632060441746 >> n=19;

>> b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 2 0]'; >> cholesky_qiu_jie(n,b) 1.0e+17 *

Columns 1 through 5 0.000000000186975

0.000000487130031

-0.000037337486618

0.000907963809187 -0.010389239660544

Columns 6 through 10 0.065117610127561

-0.236770662992548

0.490423724491871

-0.485761304739977 -0.006945225699907

Columns 11 through 15 0.375039361697087

-0.191724740660743

0.584539934693850

-1.664633268513667 0.883570199543987

Columns 16 through 19 1.944997967101147 -0.384572610629205

n=20;

b=[1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 2 0 0]'; >> cholesky_qiu_jie(n,b) 1.0e+17 * Columns 1 through 5 -0.000000005297720

0.000002201926696

-0.000135104609897

-3.195222688027543

1.831459798885819

0.003035236689788 -0.033417756031216

Columns 6 through 10 0.204643883440640

-0.727554702009794

1.433949931025033

-1.101616652841183 -1.111630914232476

Columns 11 through 15 2.722140887157841

-0.945117835365089

-0.676679844961023

-1.643697361216293 2.211735585037371

Columns 16 through 20 3.295476378249488

-7.583406568516137

5.452181213933089

-1.659242365243840 0.159333719326849

A?1利用cholesky分解求出的解的精确度高于直接，因为当n逐渐增大的过程中，

越来越接近奇异矩阵，使得计算结果的误差增大，而使用cholesky分解可以避免这种现象的产生，是计算结果更加精确。

6：（1）编制程序House(x)，其作用是对输入的向量x，输出单位向量u使得

(I?2uuT)x?x2e1。

编制Householder变换阵H?I?2uuT?Rn?n乘以A?Rn?m的程序HA，注意，你的程序并不显式的计算出H。

考虑矩阵

?1???1A???2???10?0?232224??23?e??

??37?75/2??3用你编制的程序计算H使得HA的第一列为?e1，并将HA的结果显示出来。 编制House(x)，其作用是对输入的向量x,输出单位向量U使得

(I?2uuT)x?x2e1：

House(x)通用程序： function House(x) [n,~]=size(x); e1=zeros(n,1); e1(1)=1; w=x-norm(x)*e1; U=w/sqrt((w'*w)); disp('U='); disp(U);

编制Householder变换阵H?I?2uuT?Rn?n乘以A?Rn?n的程序HA，注意，你的

程序并不显示的计算出H。

Householder变换阵H通用程序： function Householder(A) [n,~]=size(A); e1=zeros(n,1); e1(1)=1; x=A(:,1); w=x-norm(x)*e1; U=w/sqrt((w'*w)); H=eye(n)-2*(U*U'); HA=H*A; disp('HA='); disp(HA); 考虑矩阵

?1??-1A??-2???10?0?24??323?2e??

?2?37?2752??3用你编制的程序计算H使得HA的第一列为?e1的形式，并将HA的结果显示：

HA的结果显示的通用程序：

function Householder(A) [n,~]=size(A); e1=zeros(n,1); e1(1)=1; x=A(:,1);

w=x-norm(x)*e1; U=w/sqrt((w'*w)); H=eye(n)-2*(U*U'); HA=H*A; disp('HA='); disp(HA); 程序运行结果：

>> A=[1 2 3 4;-1 3 sqrt(2) sqrt(3);-2 2 exp(1) pi;-sqrt(10) 2 -3 7;0 2 7 5/2]; >> Householder(A) HA=

4.0000 -2.8311 1.4090 -6.5378 0.0000 1.3896 0.8839 -1.7805 0.0000 -1.2208 1.6576 -3.8836 0.0000 -3.0925 -4.6770 -4.1078 0 2.0000 7.0000 2.5000

7：用jacobi和Gauss-Scidel迭代求解下面的方程组，输出每一步的误差xk?x*；

?5x1?x2?3x3??2???x1?2x2?4x3?1 ??3x?4x?15x?10123?Jacobi迭代通用程序： function jacobi(A,b,x0,ep) if nargin==3 ep=1.0e-6; else if nargin< 3 error return

end end

D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\\(L+U); f=D\\b; y=B*x0+f; n=1;

while norm(y-x0)>=ep fprintf('%d步误差\\n',n) Error=y-x0; disp (Error) x0=y; y=B*x0+f; n=n+1; end

disp('近似解') disp(y)

在命令窗口输入以下命令就可以输出每一步的误差及最优近似解： >> A=[5 -1 -3;-1 2 4;-3 4 15]; >> b=[-2 1 10]'; >> x0=[0 0 0]'; >> ep=10e-8; >> jacobi(A,b,x0,ep)

近似解

-0.081967231207363 -1.803278647380326 1.131147556193974 Gauss-Scidel迭代法通用程序： function Gaussseidel(A,b,x0,ep) if nargin==3 ep=1.0e-6; else if nargin< 3 error return end end

D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=(D-L)\\U; f=(D-L)\\b; y=B*x0+f; n=1;

while norm(y-x0)>=ep fprintf('%d步误差\\n',n) Error=y-x0; disp (Error) x0=y; y=B*x0+f;

n=n+1; end

disp('x的近似最优解：') disp(y) ;

在命令窗口输入以下命令就可以输出每一步的误差及最优近似解： >> A=[5 -1 -3;-1 2 4;-3 4 15]; >> b=[-2 1 10]'; >> x0=[0 0 0]'; >> ep=10e-8;

>> Gaussseidel(A,b,x0,ep) x的近似最优解： -0.081967202756916 -1.803278585134911 1.131147515484593

8：取不同的初值用Newton迭代法以及弦截法求解x3?2x2?10x?100?0的实根，列表或者画图说明收敛性；

Newton迭代法通用程序：

%使用说明 f为符号表达式，x0为初始值 function Newton(f,x0) syms x; x=x0;

x1=x0-eval(f/diff(f)); while(norm(x1-x0)>0.000001) x=x1;

xk=x1-eval(f/diff(f)); x0=x1;

x1=xk; disp(xk); end disp('xk=') disp(xk)

弦截法通用程序：

%使用说明 f为符号表达式，x0,x1为初始值 function Xuanjiefa(f,x0,x1) syms x; x=x0; f1=eval(f); x=x1; f2=eval(f);

x2=x1-(f2/(f2-f1))*(x1-x0); while(norm(x2-x1)>0.000001) x=x1; f1=eval(f); x=x2; f2=eval(f);

xk=x2-(f2/(f2-f1))*(x2-x1); x1=x2; x2=xk; end disp('xk=') disp(xk)

Newton法：

k 0 1 2 3 xk 3 3.5102 3.4611 3.4606 f(xk) -25 2.9962 0.0301 3.1407e-06 xk?1?xk 0.5102 0.0491 5.0362e-04 k 0 1 2 3 4 xk 2 4.1333 3.5404 3.4619 3.4606 f(xk) -64 46.1178 4.8532 0.07744 2.0766e-5 xk?1?xk 2.1333 0.5929 0.0786 0.0013 k 0 1 2 3 xk 4 3.5136 3.4612 3.4606 f(xk) 36 3.1982 0.0342 4.059e-6 xk?1?xk 0.4865 0.0524 5.726e-4 弦截法：

k 0 1 2 3 4 5 xk 3 3.2 3.4879 3.4591 3.4606 3.406 f(xk) -25 -14.7520 1.6418 -0.0907 -5.1186e-04 1.6105e-07 xk?1?xk 0.2 0.2879 0.0288 0.0015 8.5667e-6 k 0 1 2 3 4 5 6 xk 2 3 3.6410 3.4428 3.4599 3.4606 3.4606 f(xk) -64 -25 11.1937 -1.0607 -0.0390 1.4420e-4 -1.9486e-8 xk?1?xk 1 0.6410 0.1983 0.0172 6.5469e-6 2.4129e-6 从上述两个表的结果我们可以对照得出，在不同的初始值的条件下运行程序都具有收敛性，但是不同的初始值有不同的收敛速度；

9：用二分法求解方程excosx?2?0在区间[0,4?]上所有根。 function er_fen_fa(f,a,b) syms x;

e=0.0000001; while((b-a)>e) m=(b+a)/2; x=a; f1=eval(f); x=m; f2=eval(f); if f1*f2>=0 a=m; else b=m; end end c=(a+m)/2; disp(c);

根据函数的特点将其定义域分为4等分分别为[0,?]，[?,2?]，[2?,3?]，[3?,4?] 分别运行程序： >> syms x

>> f=exp(x)*cos(x)+2; >> a=0; >> b=pi; >> er_fen_fa(f,a,b 近似零点： 1.8807 >> syms x

>> f=exp(x)*cos(x)+2;

a=pi; >> b=2*pi; >> er_fen_fa(f,a,b) 近似零点： 4.6941 >> syms x

>> f=exp(x)*cos(x)+2; >> a=2*pi; >> b=3*pi; >> er_fen_fa(f,a,b) 近似零点： 7.8548 >> syms x

>> f=exp(x)*cos(x)+2; >> a=3*pi; >> b=4*pi; >> er_fen_fa(f,a,b) 近似零点： 10.9955

;x2?4.6941;x2?7.8548;x4?10.9955 由上述程序可知其全部解为x1?1.880710：考虑函数f(x)?sin(?x),x?[0,1]。用等距节点作f(x)的Newton插值，求出插值多项式以及f(x)的图像，观察收敛性。

Newton插值求插值多项式的通用程序： function P=Nweton_xhazhi(X,f) syms x; p=0;

[~,n]=size(X); for i=1:n-1 F=0; w=1; sum=1; for l=1:i+1; w=w*(x-X(l)); end for j=1:i+1 x=X(j);

F=F+f(j)/(eval(diff(w))); end for t=1:i syms x;

sum=sum*(x-X(t)); end p=p+F*sum; end P=p+f(1); Disp(‘P=’);

用等距节点作f(x)的Newton插值多项式的通用程序： X=0:0.1:1; f=sin(pi*X);

O=Nweton_xhazhi(X,f); Y=simple(O); disp(Y);

plot(X,f) P=

x*(11153755408641612500000*x^9 19582893765599948750000*x^7 6457490747185549971250*x^5 452380176499959382000*x^3

+ + -

-

55768777028538456250000*x^8 256281087074564729875000*x^6 1150585140486343857745625*x^4 2327269296084898261641325*x^2

+ + + +

4715634276489610818*x - 1414847701376318005254584)

画图程序： X=0:0.1:1; x1=0:0.0001:1; f1=sin(x1*pi); syms x x0=0:0.001:1; f=sin(pi*X);

O=Nweton_xhazhi(X,f); x=x0; Y=eval(O); plot(x0,Y,'r'); hold on plot(x1,f1);

1.210.80.60.40.20-0.200.10.20.30.40.50.60.70.80.91

由上图可知：红色图形是利用插值多项式画出的，蓝色图像是利用原函数画出的，将两者对比我们可以得出结论：该插值多项式具有收敛性

11：对函数f(x)?1,x?[?5,5]，取不同的节点数n,用等距节点作Lagrange1?x2插值，观察Runge现象。

Lagrange插值通用程序：

%X为节点，y为节点对应的函数值，x0为插值节点 function Lagrange(X,y,x0) syms x; [~,n]=size(X); sum=0; for i=1:n

sum=sum+y(i)*L(X,i); end x=x0; p=eval(sum); plot(x0,p,'r-') hold on

plot(X,y,'b-');

function l=L(X,k) syms x; sum1=1; sum2=1; [~,n]=size(X); if (k==1) for i=k+1:n

sum1=sum1*(x-X(i)); sum2=sum2*(X(k)-X(i)); l=sum1/sum2; end else

for i=1:k-1

sum1=sum1*(x-X(i)); sum2=sum2*(X(k)-X(i)); end for i=k+1:n

sum1=sum1*(x-X(i)); sum2=sum2*(X(k)-X(i)); end

l=sum1/sum2; end

取不同的节点数n,用等距节点作Lagrange插值，观察Runge现象:

当n=10;

X=-5:0.1: 5;f=1./(1+X.^2); x0=-5:5:5;Lagrange(X,f,x0) x1=-5:0.0001:5; y1=1./(1+x1.^2); plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

当n=20;

X=-5:0.1:5; f=1./(1+X.^2); x0=-5:0.5:5; Lagrange(X,f,x0) x1=-5:0.0001:5; y1=1./(1+x1.^2); plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

当n=50;

X=-5:0.1:5; f=1./(1+X.^2); x0=-5:0.2:5; Lagrange(X,f,x0) x1=-5:0.0001:5; y1=1./(1+x1.^2); plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

n=100; X=-5:0.1:5; f=1./(1+X.^2); x0=-5:0.1:5;

Lagrange(X,f,x0) x1=-5:0.0001:5; y1=1./(1+x1.^2); plot(x1,y1)

1.41.210.80.60.40.20-5-4-3-2-1012345

由以上几个图形我们可以得知，插值节点越多时，插值多项式的收敛性越好； 12：令f(x)?e3xcos(?x)，考虑积分?2?0f(x)dx。区间分为50，100，200，500，

1000等，分别用复合梯形以及复合Simpson积分公式计算积分值，将数值积分的结果与精确值比较，列表说明误差的收敛性。

复合梯形公式通用程序：

%f为被积函数，n为积分区间等分的数目，a,b分别为积分下上限 function tixing(f,n,a,b) syms x; sum=0; xk=a:(b-a)/n:b; for i=2:n-1 x=xk(i);

sum=sum+eval(f); syms x; end x=a; f1=eval(f); x=b; f2=eval(f);

T=(b-a)/(2*n)*(f1+f2+2*sum); disp(T)

复合sipson公式： function Simpson(f,n,a,b) syms x; sum=0; sum1=0; xk=a:(b-a)/n:b; for i=1:n-1

x=(xk(i)+xk(i+1))/2; sum=sum+eval(f); syms x; end for i=2:n-1 x=xk(i);

sum1=sum1+eval(f); syms x; end x=a;

f1=eval(f); x=b; f2=eval(f);

T=(b-a)/(6*n)*(f1+f2+2*sum1+4*sum); disp(T)

利用matlab内置函数int()求出该积分式的真是值为： >> syms x

>> int(exp(3*x)*cos(pi*x),0,2*pi) ans =

3.5232e+07

利用符合simpson公式计算如下表：

n 50 100 200 500 1000 2000 3000 5000 近似值FK 2.3147e+07 2.9066e+07 3.2156e+07 3.4010e+07 3.4623e+07 3.4928e+07 3.5030e+07 3.5111e+07 误差F-Fk 1.2085e+07 6.1664e+06 3.0761e+06 1.2228e+06 6.0960e+05 3.0430e+05 2.0275e+05 1.2159e+05

利用复合梯形公式计算如下表：

n 50 100 200 500 1000 2000 3000 5000 近似值FK 2.3478e+07 2.9054e+07 3.2139e+07 3.4005e+07 3.4622e+07 3.4928e+07 3.5030e+07 3.5111e+07 误差F-Fk 1.1755e+07 6.1787e+06 3.0939e+06 1.2273e+06 6.1085e+05 3.0463e+05 2.0290e+05 1.2165e+05

由上述两个表格的值我们可以得出，其误差随着n的增大近似值与真实值之间的误差逐渐减少，说明了误差具有收敛性。

13：分别用2点,3点以及5点的Gauss型积分公式计算如下定积分：

(1)?1x21?x2?1dx(2)?2?0sinxdx x解：利用Gauss型求积公式得出：

n 0 1 xk 0 ?0.5773502692 2 ?0.77459666920 高斯型求积公式部分节点、系数表 n Ak xk 2 5 ?0.9324695142 1 ?0.6612093865 ?0.23861918616 ?0.94910791230.5555555556 0.8888888889?0.74153118563 4 ?0.8611363116 ?0.3399810436?0.90617984590.3478548451 0.65214515290.2369268851 7 ?0.5384693101 0 0.4786286705 0.5688888890.40584515140?0.9602898566?0.7966664774 ?0.5255324099?0.1834346425 Ak 0.17132449240.3607615730 0.46791393460.12948496620.2797053915 0.38183005050.41795918370.10122853630.2223810345 0.31370664590.36268378342点Gauss积分公式求解时，n=1;查找上表得x0,x1,A0,A1，把他们都代入，可得

?1x21?x2?1dx?(0.5773502692)21?(0.5773502692)2?(?0.5773502692)21?(?0.5773502692)2?0.8165

3点Gauss积分公式求解时，n=2；查找上表可得x0,x1,x2,A0,A1,A2，把他们都代入，可得：

?1x21?x2?1dx?0.5555555556?(0.7745966692)21?(0.7745966692)222?

0.5555555556?(?0.7745966692)1?(?0.7745966692)?1.05415点Gauss积分公式求解时，n=4;查找上表可得x0,x1,x2,x3,x4,A0,A1,A2,A3,A4，把他们都代入，可得：

?1x22?11?x0.53846931012dx?2?0.2369268851??1.24950.906179845921-0.90617984592?2?0.4786286705?1-0.5384693101

2

?2?02?sinx1sin(?(t?1))sinx2??02??0dx，令x?t???t??则?dx??dt，于是利用

0-1x22xt?1上表计算Gauss积分

2点Gauss积分公式求解时，n=1;查找上表得x0,x1,A0,A1，把他们都代入，可得

1sin(?(t?1))sinxsin(?(1?0.5773502692))dx?dt???0x?-1t?11?0.5773502692

sin(?(1?0.5773502692))?1.68121?0.57735026922?3点Gauss积分公式求解时，n=2；查找上表可得x0,x1,x2,A0,A1,A2，把他们都代入，可得：

1sin(?(t?1))sinxsin(?(1?0.7745966692))dx?dt?0.5555555556??0x?-1t?11?0.7745966692

sin(?(1?0.7745966692))?0.5555555556??1.39951?0..77459666922?5点Gauss积分公式求解时，n=4;查找上表可得x0,x1,x2,x3,x4,A0,A1,A2,A3,A4，把他们都代入，可得：

1sin(?(t?1))sinx51??0xdx??-1t?1dt?0.2369268859))sin(?(1?0.9061798459))??sin(?(1?0.90617984????1?0.90617984591?0.9061798459??2?01))sin(?(1-0.5384693101))??sin(?(1?0.538469310.4786286705????

1?0.53846931011-0.5384693101???1.418114：考虑微分方程初值问题：

1?dx?(tx?x2)?2 ?dt(t?1)?x(0)?2?分别用Euler法，改进的Euler法，Runge-Kutta法求解该方程。分别去步长为0.1，0.01，0.001，计算到x(1),画图说明结果。

Euler法：

本文来源：https://www.bwwdw.com/article/vne3.html