西南交大MATLAB上机题解.docx

西南交大MATLAB上机题解.docx

《西南交大MATLAB上机题解.docx》由会员分享，可在线阅读，更多相关《西南交大MATLAB上机题解.docx（19页珍藏版）》请在冰豆网上搜索。

西南交大MATLAB上机题解

序言

根据实习内容及自身学习情况，文章采用的计算软件是MATLAB。

通过在该软件M文件中编写相应程序来解决实习内容中的问题。

MATLAB软件的发展已经30多年，过程中不断的版本升级和功能完善，使得该软件具有非常强大的编程和计算能力，也使其成为当前世界上应用最普遍的计算软件之一。

作为一款不断发展的计算软件，MATLAB具有6个非常突出的特点：

（1）友好的工作平台和编程环境；

（2）简单易用的编程语言；（3）强大的科学计算数据处理能力；（4）出色的图形处理功能；（5）应用广泛的模块集合工具箱；（6）实用的程序接口和发布平台。

这些特点保证了MATLAB软件能够长期在工程界和学术界得到应用。

科学计算长期以来都是MATLAB软件擅长的部分。

MATLAB的科学计算环境不仅提供了大量用于计算的函数库、工具箱，而且还提供了功能强大的计算绘图能力，其主要的优点有以下几个方面：

（1）MATLAB几乎集成了当前所有科学研究领域和工程计算领域的算法，使得软件的科学计算能力非常的强大，同时该软件默认的计算精度是双精度，保证了计算结果的高准确性。

（2）MATLAB将高性能的数值计算和可视化集成在一起，同时提高大量的内置函数和开放的程序和数据接口，使得该软件能够在科学计算、控制系统与信息处理等领域的分析、仿真和设计大展身手。

（3）MATLAB包含各种能够进行常规运算的工具箱，同时具有编程计算功能。

（4）MATLAB能够将计算的结果绘制成二维和三维图形，增强了计算结果的可视化程度。

（5）MATLAB的编程语言是一种解释执行的脚本语言，简单易学。

文章采用MATLAB软件完成数值分析的上机实习时基于以下几个方面的考虑：

（1）前期对MATLAB软件有一定程度的了解，能够利用该软件做一些简单的矩阵计算。

（2）MATLAB的编程语言和C++的编程语言有很多相通的地方，而本科时候学过C++编程，对一些基础的编程掌握比较好，容易实现MATLAB编程计算。

（3）MATLAB软件应用非常的广泛，作为一个工科生，非常有必要掌握该软件的一些使用功能。

第一题

写出对一般的线性方程组通用的Gauss消元,Gauss-Seidel迭代程序。

并以下面的线性方程组为例进行计算，讨论所得到的计算结果是否与理论一致。

（1）

（2）

（3）

1.Gauss消元法：

此法是将待求解的线性代数方程组中各个参数表中成对应的矩阵和向量，对矩阵进行有目的的消元，使原来的矩阵经过消元变成一个上三角矩阵，以此求得方程组的最后一个方程的解，然后从后往前回代，逐步求得各个未知数的解。

一般线性方程组通用的Gauss消元法的迭代程序见附录一。

题例中的求解方式是首先在MATLAB软件中编写M文件，保存运行。

然后在MATLAB的CommandWindow中输入对应的矩阵和向量，进行求解。

题目1在CommandWindow中的程序：

a=[6,2,-1;1,4,-2;-3,1,4];

b=[-3;2;4];

x=Gauss（a,b）

运行结果：

x=[;;]

题目2在CommandWindow中的程序：

a=[1,,;,1,;,,1];

b=[3;2;1];

x=Gauss（a,b）

运行结果：

x=[;;]

题目3在CommandWindow中的程序：

a=[1,3;-7,1];

b=[4;6];

x=Gauss（a,b）

运行结果：

x=[;]

2.Gauss-Seidel迭代法

Gauss-Seidel迭代算法类似Jacobi迭代算法，Jacobi迭代算法的一般迭代格式为：

式

（1）

在Jacobi迭代算法收敛时，式

（1）中第一个方程计算得到的

马上投入到第二方程中使用，把前两个方程已经得到的

，

马上投入到第三个方程中使用，如此循环，通过这种计算线性方程组解的迭代格式称为Gauss-Seidel迭代，其迭代格式为：

式

（2）

矩阵形式如式（3）所示：

式（3）

解得

，得到的结果如式（4）所示：

式（4）

一般线性方程组通用的Gauss-Seidel迭代法的迭代程序见附录二。

题例中的求解方式是首先在MATLAB软件中编写M文件，保存运行。

然后在MATLAB的CommandWindow中输入对应的矩阵和向量，进行求解。

由于Gauss-Seidel迭代法计算时要首先考虑迭代收敛与否的问题，所以在M文件中，通过语言控制，首先对求解的线性方程组的矩阵进行相关计算，得到迭代阵B的谱半径值，如果谱半径值R<1，则进行Gauss-Seidel迭代，否则不进行迭代。

题目1在CommandWindow中的程序：

A=[6,2,-1;1,4,-2;-3,1,4];

b=[-3;2;4];

x0=[0;0;0];

x=GaussSeidel（A,b,x0）

经过计算，得到题目1的迭代阵B的谱半径R=<1，所以能够得到收敛解。

经过41次迭代后，得到该线性方程组的解为：

x=[;;]

题目2在CommandWindow中的程序：

A=[1,,;,1,;,,1];

b=[3;2;1];

x0=[0;0;0];

x=GaussSeidel（A,b,x0）

经过计算，得到题目2的迭代阵B的谱半径R=<1，所以能够得到收敛解。

经过12次迭代后，得到该线性方程组的解为：

x=[;;]

题目3在CommandWindow中的程序：

A=[1,3;-7,1];

b=[4;6];

x0=[0;0];

x=GaussSeidel（A,b,x0）

经过计算，得到题目3的迭代阵B的谱半径R=21>1，所以该线性方程组不能得到用Gauss-Seidel迭代方法求得收敛解。

通过使用上面三个例题求解的过程发现，在MATLAB中编写的M文件能够很好地求解线性方程组，Gauss消元法能够得到上面三个例题的全部收敛解，而Gauss-Seidel虽然无法求解到题目3的收敛解，但是整个M文件程序的使用性能还是非常的不错，能够真实反映出线性方程组系数矩阵内在的一些性质。

第二题

问题：

给定初值问题

（精确解为

），用Runge-Kutta4阶算法按步长

求解，分析其中遇到的现象及问题。

龙格-库塔法（Runge-Kutta）经常用来解决常微分方程数值问题，与欧拉（Euler）法、线性多步法、Adams内插公式等其他几种常微分方程数值求解方法在应用上各有千秋。

龙格-库塔法的推导是基于泰勒（Taylor）展开方法，故而此法要求求解的微分方程解具有较好的光滑性质。

根据推导得到四阶龙格-库塔法的具体格式如式（5）所示：

其中：

式（5）

求解本题微分方程的程序保存在M文件中，具体的内容见附录三。

当步长h=时，运行程序能够得到相应的结果，将该结果绘制成图形，如图

（1）所示：

图

（1）

当h=时，结果对应的图形如图

（2）所示：

图

（2）

当h=时，结果对应的图形如图（3）所示：

图（3）

第三题

问题：

给定函数

，及节点

，利用高次多项式、分段线性多项式进行插值，画图比较几种插值多项式与

的图形，并说明其中的数学现象。

利用Newton插值法对该函数进行插值，通用的Newton插值法程序保存在M文件中，具体程序清单见附录四。

实际插值时，先将对应的M文件打开，在MATLAB的CommandWindow中输入下面的程序：

x=[-5:

1:

5];

y=1./（1+x.^2）;

x0=[-5:

:

5];

y0=Lagrange（x,y,x0）;

y1=1./（1+x0.^2）;

gridon;

plot（x0,y0,'--r'）;

holdon;

plot（x0,y1,'--b'）;

xlabel（'自变量--x'）;

ylabel（'因变量--y'）;

legend（'Lagrange插值图形','函数准确图形'）;

title（'Lagrange插值图形和函数准确图形'）;

gridon

最终得到如图（4）的图形：

图（4）

利用分段线性插值多项式进行插值计算，程序保存在M文件中，具体程序清单见附录四。

实际插值时，先将对应的M文件打开，在MATLAB的CommandWindow中输入下面的程序：

x=[-5:

:

5];

y=1./（1+x.^2）;

u=[-5:

1:

5];

v=pline（x,y,u）;

gridon;

plot（x,y,'--r'）;

holdon;

plot（u,v,'--b'）;

xlabel（'自变量--x'）;

ylabel（'因变量--y'）;

legend（'pline插值图形','函数准确图形'）;

title（'pline插值图形和函数准确图形'）;

gridon

最终得到如图（5）的图形：

图（5）

理论上而言，对于一个给定区间[a,b]，根据节点分布采用插值多项式得到函数的近似值时，插值多项式的次数越高，逼近

值的精度就越高。

Newton法插值得到的多项式

次数明显高于分段线性插值，所以Newton插值法逼近

值的精度要高于分段线性插值法，但是通过上面两个实验获得的图形可以发现，事实并非如此，图（4）和图（5）展示出了Newton插值法在区间的首位处逼近

值的精度明显低于分段线性插值，从而得出高次多项式插值的精度并不是在区间上任何地方都高于分段线性的结论。

造成本题出现该问题产生的原因是：

函数

在区间[-5,5]上的各阶导数存在，但是由n个节点所构成的Lagrange插值多项式在全区间内并非都收敛。

第四题

问题：

利用牛顿法求方程

于区间

的根，考虑不同初值下牛顿法的收敛情况。

迭代法是求解非线性方程（组）的主要方法，而Newton法是最有效的迭代方法之一，它在单根附近具有较高阶的收敛速度，统统是还能够求解到方程的复根（此时初值需要时复数）。

通常是在方程

的根

附近，用

的某个近似值去代替

，再令这个近似式为零，解出的根就作为根

的新的近似值。

基于这样的一种思想，可以推导得到Newton法的计算公式，首先设

的近似根为

，函数

可以在

附近作Taylor展开，得到如式（6）：

式（6）

利用式（6）的线性部分作为函数

的近似，得到如式（7）的公式：

式（7）

求解式（7），得到根值为

，即（8）：

式（8）

由此推断得到Newton迭代公式为式（9）：

式（9）

针对此题的MATLAB程序见附录六，由于题中给定根值区间为

，故首先取初值

在MATLAB中CommandWindow窗口中输入以下应用程序：

f=inline（'x-log（x）-2'）;

df=inline（'1-1./x'）;

Newton（f,df,,

运行得到结果为：

第0次迭代（初值）x[0]=

第1次迭代x[1]=

第2次迭代x[2]=

第3次迭代x[3]=

第4次迭代x[4]=

ans=

改变初值

，运行得到结果为：

第0次迭代（初值）x[0]=

第1次迭代x[1]=

第2次迭代x[2]=

第3次迭代x[3]=

ans=

改变初值

，运行得到结果为：

第0次迭代（初值）x[0]=

第1次迭代x[1]=

第2次迭代x[2]=

第3次迭代x[3]=

第4次迭代x[4]=

ans=

改变初值

，运行得到结果为：

第0次迭代（初值）x[0]=

第1次迭代x[1]=

第2次迭代x[2]=

第3次迭代x[3]=

第4次迭代x[4]=

第5次迭代x[5]=

ans=

改变初值

，运行得到结果为：

第0次迭代（初值）x[0]=

第1次迭代x[1]=

第2次迭代x[2]=

第3次迭代x[3]=

第4次迭代x[4]=

第5次迭代x[5]=

ans=

根据不同初值的迭代过程可以发现，不同的初值对应的迭代次数并不完全相同，距离最终根的距离越近，迭代的次数相对越少。

经过几次取值，计算结果都没有出现根不收敛的情况，所以无法分析初值对根最终收敛的影响。

附录

第一题程序清单

Gauss迭代的M文件源程序

function[x]=Gauss（a,b）

n=length（a）;

x=zeros（n,1）;

a=[ab];

fork=1:

n-1

max=k;

fori=k+1:

n

ifa（i,k）>a（max,k）

max=i;

end

end

temp=a（k,k:

n+1）;

a（k,k:

n+1）=a（max,k:

n+1）;

a（max,k:

n+1）=temp;

fori=k+1:

n

a（i,k）=-a（i,k）/a（k,k）;

a（i,k+1:

n+1）=a（i,k+1:

n+1）+a（i,k）*a（k,k+1:

n+1）;

end

end

x（n,1）=a（n,n+1）/a（n,n）;

fori=n-1:

-1:

1

sum=0;

forj=i+1:

n

sum=sum+x（j,1）*a（i,j）;

end

x（i,1）=（a（i,n+1）-sum）/a（i,i）;

end

Gauss-Seidel迭代的M文件源程序：

functionx=GaussSeidel（A,b,x0）

D=diag（diag（A））;

L=-tril（A,-1）;

U=-triu（A,1）;

B=（D-L）\U;

f=（D-L）\b;

x=zeros（length（x0）,1）;

i=0

x0=x;

x=B*x+f;

R=max（abs（eig（B）））

ifR<1

fori=1:

100

ifx0==x;

i=i+1

break;

else

i=i+1;

x0=x;

x=B*x0+f;

end

end

else

fprintf（'不能够用Gauss-Seidel迭代法计算’）

end

end

第二题程序清单

Runge-Kutta4阶算法程序清单

clear;

F='-20*x';

a=0;b=2;

h=;

n=（b-a）/;

X=a:

:

b;

Y=zeros（1,n+1）;

Y

（1）=1;

fori=1:

n

x=X（i）;

y=Y（i）;

K1=eval（F）;

x=x+h/2;

y=y+h*K1/2;

K2=eval（F）;

y=Y（i）+h*K2/2;

K3=eval（F）;

x=X（i）+h;

y=Y（i）+h*K3;

K4=eval（F）;

Y（i+1）=Y（i）+h*（K1+2*K2+2*K3+K4）/6;

end

%准确值

temp=[];

f=dsolve（'Dy=-20*x','y（0）=1','x'）;

df=zeros（1,n+1）;

fori=1:

n+1

temp=subs（f,'x',X（i））;

df（i）=double（vpa（temp））;

end

disp（'步长四阶经典R-K法准确值'）;

disp（[X',Y',df']）;

figure;

plot（X,df,'k*',X,Y,'--r'）;

gridon;

title（'四阶经典R-K法解常微分方程'）;

legend（'准确值','四阶经典R-K法'）;

第三题程序清单

高次多项式插值法（Newton法）程序清单

functionyy=Lagrange（x,y,xi）

m=length（x）;

n=length（y）;

ifm~=n,error（'向量x和y的长度必须一致'）;

end

s=0;

fori=1:

n

z=ones（1,length（xi））;

forj=1:

n

ifj~=i

z=z.*（xi-x（j））/（x（i）-x（j））;

end

end

s=s+z*y（i）;

end

yy=s;

MATLAB中CommandWindow输入的应用程序清单

x=[-5:

1:

5];

y=1./（1+x.^2）;

x0=[-5:

:

5];

y0=Lagrange（x,y,x0）;

y1=1./（1+x0.^2）;

gridon;

plot（x0,y0,'--r'）;

holdon;

plot（x0,y1,'--b'）;

xlabel（'自变量--x'）;

ylabel（'因变量--y'）;

legend（'Lagrange插值图形','函数准确图形'）;

title（'Lagrange插值图形和函数准确图形'）;

gridon

分段线性多项式插值法程序清单

functionv=pline（x,y,u）

delta=diff（y）./diff（x）;

n=length（x）;

k=ones（size（u））;

forj=2:

n-1

k（x（j）<=u）=j;

end

s=u-x（k）;

v=y（k）+s.*delta（k）;

MATLAB中CommandWindow输入的应用程序清单

x=[-5:

:

5];

y=1./（1+x.^2）;

u=[-5:

1:

5];

v=pline（x,y,u）;

gridon;

plot（x,y,'--r'）;

holdon;

plot（u,v,'--b'）;

xlabel（'自变量--x'）;

ylabel（'因变量--y'）;

legend（'pline插值图形','函数准确图形'）;

title（'pline插值图形和函数准确图形'）;

gridon

第四题程序清单

Newton法求解非线性方程组程序清单

functionx=Newton（f,df,x0,e,N）

ifnargin<5

N=500;

end

ifnargin<4

e=1e-4;

end

x=x0;

x0=x+2*e;

k=0;

fprintf（'第%2d次迭代（初值）x[%2d]=%\n',k,k,x）

whileabs（x0-x）>e&k

k=k+1;

x0=x;

x=x0-feval（f,x0）/feval（df,x0）;

fprintf（'第%2d次迭代x[%2d]=%\n',k,k,x）

end

ifk==N,fprintf（'迭代次数已经达到上限值'）;

end;