课程实验及解决方案

?数值计算?实验指导书

要求：

（1） 用Matlab语言或你熟悉的其它算法语言语言编程序，使之

尽量具有通用性．

（2） 上机前兖分准备，复习有关算法，写出计算步骤，反复查对

程序．

（3） 完成计算后写出实验报告报告，内容包括：计算机型号及所

用算法语言，CPU时间，结果分析和小结等． （4） 熟练掌握Marlab的使用．

实验－

实验目的：掌握Matlab在数值计算求解中的基本方法． 实验要求：会使用Matlab工具．并能用它求解以下问题问题． 1．

分别用库函数quad8和quad求以下定积分：

??sinxdx

12x2． 库函数eig求下列矩阵的特征值和特征向量：

?4?11?? A???13?2????1?23??3． 作下列函数在指定区间上的图形，并用库函数roots求它的

根：

f(x)?x?4x?6x?1 x?[0,4]

32

1

实验二

实验目的：掌握三次样条插值的三弯矩法． 实验要求：运用三弯矩法编制程序求下列问题的解． 实验题： 设已知数据如下：

x i0.2 0. 4 0.6 0.8 1.0 f(xi) 0.9798652 0.9177710 0.8080348 0.6386093 0.3843735 求f(x)的三次样条插值函数 s(x)，满足： （1） 自然边界条件

s''(0.2)?s(1.0)?0;

''（2） 第一种边界条件

s(0.2)?0.20271,s(1.0)?1.55741.要求输出

''用追赶法解出的弯矩向量

(M,M01,?,M4)和

s(0.2?0.1i)(i?0,1,?,8)的值，并画出y?s(x)的图形．

实验三

实验目的：复化求积公式及高斯－勒让德公式计算定积分. 实验题：

（1）ln2??2?2实验要求：

（1） 分别用复化梯形公式，复化Simpson公式求解，要求绝对误

?7差限为???10．

311 (2) dx??4dx ?2201?xx?1112（2） 用三点的五点的Gauss-Legendre求积公式求解上述积分，

2

并输出I的近似值． （3） 分析比较各种计算结果．

实验四

实验目的：迭代法求解非线性方程的根． 实验题： 实验要求：

（1） 用你自己设计的一种线性收敛迭代法求方程的根，然后

用斯蒂芬森加速迭代重新计算．

（2） 用牛顿法求方程的根，输出迭代初值．各次迭代值及迭

代次数，并与（1）比较．

x3?3x?2???0

x实验五

实验目的：观察的理解高斯消元过程中出现小主元即akk很小时，引

起方程组解数值不稳定性．

实验题：求解线性方程组AX?b．

?15?0.3?1059.143?5.291?6.130?1 A???11.295?121??(k)?59.17?1???46.78?2?? ，b???1?2????21????实验要求：

（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的． （2） 用高斯列主元消去法求得L和U及解向量x1,x2?R．

3

4（3） 用不选主元的高斯消去法求得L和U及解向量x1,x2?R． （4） 观察小主元并分析对计算结果的影响．

4实验六

实验目的：认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性

质对收敛速度的影响．

实验题：用迭代法求解方程组AX?b,其中A?R20?20，它的每条对角线

元素是常数．

??3?1???2??1?A??4????????1212141?2?3??14??????? 1???4?1??2??3??3?141?2??121?4???3?12?实验要求：

（1） 选取不同的初始向量x和不同的方程组右端向量b，给定

迭代误差要求．用Jacobi和Gauss-Seidel迭代法计算，观测得到的迭代向量序列是否收敛，若收敛，记录迭代次数，分析计算结果并得出你的结论．

（2） 取定右端向量b和初始向量x，将的主对角线元素成倍增

长若干次，非主对角线元素不变，每次用Jacobi迭代法计算，要求迭代误差满足

(0)(0)

x(k?1)?x(k)??10，比较收敛速度，

?54

分析现象并得出你的结论．

实验七

实验目的：求初值问题的常微分方程的数值解 实验题：

给定初值问题为

y?'1??,1?x?2,?y2 ?xx?y(1)?1.?实验要求：

（1） 用改进欧拉法（取h=0.05）求它的数值解，并输出

x?0.1i(i?0,1,2,?,10)的数值解y

ii（2） 用四阶R-K方法（取h=0.1）求它的数值解，并输出

x?0.1i(i?0,1,2,?,10)的数值解y

ii（3）

分析结果，对两种算法进行比较．

综合实验

实验目的：由实际问题建立数学模型，然后会用相应的数值算法求解． 实验题：生态环境学家在研究自然界中两个生物种群数目变化时得到一组常微分方程。

假设有两种生物（例如一种是蓝鲸，另一种是南极磷虾），前者在时刻t时的数量为x1(t)，后者在时刻t时的数目为x2(t)，并假设它们都是t的连续可微函数。蓝鲸是以磷虾为主要食物的。当没有食物来

5

td?xxb?2xa?2xd??td???2x1xd?1xc??1xd?21源时蓝鲸数目会减少，其减少速度与当时蓝鲸的数目成线性关系，即

dx1??cx1(t) （1） . dt当有食物来源时，蓝鲸的数目会增加。增加的速度和它捕食的数目有关，即

dx1= dx1(t) x2(t) （2） . dt合并（1）和（2），得到蓝鲸变化速度满足的微分方程

dx1??cx1(t)? dx1(t) x2(t). (3) dt同样，在没有蓝鲸时，磷虾的增加速度满足

dx2=ax(t). (4) dt考虑到被捕食情况，则磷虾的数目满足

dx2=ax(t)－bx1(t) x2(t). (5) dt合并（3）和（5），得到著名的Lotka-Volterra方程

（6）

其中a,b,c,d均为正常数。

（6）是一个非线性常微分方程组，不可能有解析解。

实验要求：假设a?1.2,b?0.6,c?0.8,d?0.3,而且初始值为x1(0)=2,

x2(0)=1.

1）

用不同的数值方法求解（6）。把x1(t) 和x2(t)画在同一张

图上并给以解释。

6

2） 把（6）的两个方程相除，得到

dx1?cx1?dx1x2 （7） ?dx2ax2?bx1x2用不同的数值方法解出x2~x1之间的函数关系。并把它画在以x1,x2为坐标的图上，对所得结果加以解释。

附1：部分数值分析常用实验程序

1）Guass列主元消元法解线性方程组：

//高斯列主元消去法求解非线性方程组源程序清单 # define NUMBER 20 # include \

float A[NUMBER][NUMBER+1] ,ark; int flag,n;

exchange(int r,int k) { int i;

for(i=1;i<=n+1;i++) A[0][i]=A[r][i];

for(i=1;i<=n+1;i++) A[r][i]=A[k][i];

for(i=1;i<=n+1;i++) A[k][i]=A[0][i];} float max(int k) { int i;

float temp=0;

for(i=k;i<=n;i++)

if(fabs(A[i][k])>temp) { temp=fabs(A[i][k]); flag=i;} return temp; } main()

{ float x[NUMBER],b[NUMBER]; int r,k,i,j;

me: printf(\用Gauss列主元消元法解线性方程组\ printf(\ scanf(\

printf(\现在输入系数矩阵A和向量b:\\n\\n\

7

for(i=1;i<=n;i++) for(j=1;j<=n+1;j++)

{ if(j==n+1) printf(\请输入b%d:\ else printf(\请输入a%d%d:\ scanf(\ } k=1; label:

/****************************************************/ ark=max(k);

/*if(fabs(A[1][k])>=fabs(A[2][k])) { ark=fabs(A[1][k]);flag=1;} else { ark=fabs(A[2][k]);flag=2;} if(fabs(A[3][k])>ark)

{ ark=fabs(A[3][k]);flag=3; } */

/**************************************************/ if(ark==0)

{ printf(\不是合法的方程组!\ else if(flag!=k) exchange(flag,k); for(i=k+1;i<=n;i++) for(j=k+1;j<=n+1;j++)

A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k]; if(k

{ k++; goto label;} else

{ x[n]=A[n][n+1]/A[n][n]; for( k=n-1;k>=1;k--) {float me=0;

for(j=k+1;j<=n;j++) { me=me+A[k][j]*x[j];}

x[k]=(A[k][n+1]-me)/A[k][k]; } }

for(i=1;i<=n;i++)

printf(\ goto me; }

/////////////////////////////////////////////////////////////////// 2）//插值源程序

# define f(x) 1/(1+25*x*x)

# define f1(x) -50*x/((1+25*x*x)*(1+25*x*x)) # include \float xi[11],M[11];

8

int N=10; float h=0.2;

float L(int i,float x) { float ll=1.0;int j; for( j=0;j<=N;j++) { if(i==j) continue; else ll=ll*(x-xi[j]); }

for(j=0;j<=N;j++) { if(i==j) continue;

else ll=ll/(xi[i]-xi[j]); }

return ll; }

float S(int N,float x) { float ss=0; int j; for( j=0;j<=N;++j) ss=ss+f(xi[j])*L(j,x); return ss; }

float S1(float x0,float x1,float x) { float ss;

ss=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1); return ss;}

float chashan(float x0,float x1,float x2)

{ return ((f(x2)-f(x0))/(x2-x0)-(f(x1)-f(x0))/(x1-x0))/(x2-x1); }

float S3( float xj,float xjp1 ,float Mj,float Mjp1,float x) { float ss1,ss2,ss3,ss4,ss;

ss1=(xjp1-x)*(xjp1-x)*(xjp1-x)/(6*h)*Mj; ss2=(x-xj)*(x-xj)*(x-xj)/(6*h)*Mjp1; ss3=(x-xj)/h*(f(xjp1)-h*h/6*Mjp1); ss4=(xjp1-x)/h*(f(xj)-h*h/6*Mj); ss=ss1+ss2+ss3+ss4; return ss; }

main()

{ int i,j;int control;float help1,help2,helpe,helpb; float d[11]; int cchar=0; float x0,xo1[10],xo2[10]; x0=0;

for( i=0;i<=N;i++) xi[i]=-1.0+i*2.0/N; for( j=0;j<=9;++j)

9

{ xo1[j]=0.06+0.1*j; xo2[j]=-0.06-0.1*j; } label:

printf(\现象的克服\printf(\多项式插值\printf(\分段线性插值\

printf(\三次样条函数插值\printf(\退出\printf(\请选择:\ scanf(\switch(control)

{case 1: //下面是多项式插值 for(i=0;i<=9;i++){

printf(\ printf(\ printf(\ printf(\ help1=f(xo1[i]);help2=S(N,xo1[i]); helpe=help1-help2;

printf(\helpb=helpe/help1*100.0;

printf(\相对误差=%9.6f%% \ } break;

case 2: // 下面是分段线性插值 for(j=0;j<=9;j++) for(i=0;i<=N;i++)

{if(xo1[j]>xi[i] && xo1[j]

help2=S1(xi[i],xi[i+1],xo1[j]); helpe=help1-help2;

helpb=helpe/help1*100.0;

printf(\ printf(\ printf(\printf(\相对误差=%9.6f%% \ }

if(xo2[j]>xi[i] && xo2[j]

help2=S1(xi[i],xi[i+1],xo2[j]);

printf(\ printf(\ }

10

}

break;

case 3: //以下是三次样条函数插值

d[0]=6/h*((f(xi[1])-f(xi[0]))/h-f1(xi[0])); printf(\ for(j=1;j<=N-1;j++)

{ d[j]=6*chashan(xi[j-1],xi[j],xi[j+1]); printf(\ }

d[N]=6/h*(f1(xi[N])-(f(xi[N])-f(xi[N-1]))/h); printf(\

printf(\请用Gauss列主元法解十一阶三对角方程求出M0至M10 !\

printf(\然后把M0至M10输入! 准备好后按1:\scanf(\ if(cchar==1)

{ for(j=0;j<=10;j++)

{ printf(\请输入M%d:\ scanf(\ } }

for(j=0;j<=9;j++) for(i=0;i<=N;i++)

{if(xo1[j]>xi[i] && xo1[j]

help2=S3(xi[i],xi[i+1],M[i],M[i+1],xo1[j]); helpe=help1-help2;

helpb=helpe/help1*100.0;

printf(\ printf(\ printf(\printf(\相对误差=%9.6f%% \ }

if(xo2[j]>xi[i] && xo2[j]

help2=S3(xi[i],xi[i+1],M[i],M[i+1],xo2[j]); printf(\ printf(\ } }

11

附2：综合实验参考 （1） 设计思路：

根据欧拉方法、改进的欧拉方法和经典的4级4阶Runge-Kutta法来求解。 （2） 函数文件： 公用的三个函数文件： fun.m函数文件： function y=fun(x1,x2) if -0.8*x1+0.3*x1*x2~=0

y=(1.2*x2-0.6*x1*x2)/(-0.8*x1+0.3*x1*x2); end

fun1.m函数文件： function y=fun1(x) y=-0.8*x(1)+0.3*x(1)*x(2); fun2.m函数文件： function y=fun2(x) y=1.2*x(2)-0.6*x(1)*x(2); 第一问的三种方法求解文件： euler方法： function euler(n,h) x(1)=2; x(2)=1;

12

t(1)=0.0000; p(1)=x(1); q(1)=x(2); for i=1:n

f(1)=fun1(x); f(2)=fun2(x); y=x+h*f; x=y; t(i+1)=i*h; p(i+1)=x(1); q(i+1)=x(2); end [t',p',q'];

plot(t,p,'r',t,q,'k'); xlabel('时刻t'); ylabel('数量x'); text(0,2,'蓝鲸x1'); text(0,1,'磷虾x2'); 改进的euler方法： function impeuler(n,h) x(1)=2; x(2)=1;

13

t(1)=0.0000; p(1)=x(1); q(1)=x(2); for i=1:n

f(1)=fun1(x); f(2)=fun2(x); g(1)=fun1(x+h*f); g(2)=fun2(x+h*f); x=x+1/2*h*(f+g); t(i+1)=i*h; p(i+1)=x(1); q(i+1)=x(2); end [t',p',q'];

plot(t,p,'r',t,q,'k'); xlabel('时刻t'); ylabel('数量x'); text(0,2,'蓝鲸x1'); text(0,1,'磷虾x2');

经典的4级4阶Runge-Kutta法： function rk4(n,h) x(1)=2;

14

x(2)=1; t(1)=0.0000; p(1)=x(1); q(1)=x(2); for i=1:n

s1(1)=h*fun1(x); s1(2)=h*fun2(x); s2(1)=h*fun1(x+1/2*s1); s2(2)=h*fun2(x+1/2*s1); s3(1)=h*fun1(x+1/2*s2); s3(2)=h*fun2(x+1/2*s2); s4(1)=h*fun1(x+s3); s4(2)=h*fun2(x+s3); x=x+1/6*(s1+2*s2+2*s3+s4); t(i+1)=i*h; p(i+1)=x(1); q(i+1)=x(2); end [t',p',q'];

plot(t,p,'r',t,q,'k'); xlabel('时刻t'); ylabel('数量x');

15

text(0,2,'蓝鲸x1'); text(0,1,'磷虾x2'); 第二问的三种方法求解文件： euler方法：

function deuler(n,h) x1=2:h:n;m=max(size(x1)); x2(1)=1; for i=2:m

x2(i)=x2(i-1)+h*fun(x1(i-1),x2(i-1)); end

plot(x1,x2,'-r');

title('蓝鲸和磷虾的数量关系'); xlabel('蓝鲸的数量'); ylabel('磷虾的数量'); 改进的euler方法： function dimpeuler(n,h) x1=2:h:n;m=max(size(x1)); x2(1)=1; for i=2:m

k1=fun(x1(i-1),x2(i-1)); k2=fun(x1(i),x2(i-1)+h*k1); x2(i)=x2(i-1)+h*(k1+k2)/2;

16

end

plot(x1,x2,'-r');

title('蓝鲸和磷虾的数量关系'); xlabel('蓝鲸的数量'); ylabel('磷虾的数量');

经典的4级4阶Runge-Kutta法： function drk4(n,h) x1=2:h:n;m=max(size(x1)); x2(1)=1; for i=2:m

k1=fun(x1(i-1),x2(i-1));

k2=fun(x1(i-1)+h/2,x2(i-1)+h*k1/2); k3=fun(x1(i-1)+h/2,x2(i-1)+h*k2/2); k4=fun(x1(i-1)+h,x2(i-1)+h*k3); x2(i)=x2(i-1)+h*(k1+2*k2+2*k3+k4)/6; end

plot(x1,x2,'-r');

title('蓝鲸和磷虾的数量关系'); xlabel('蓝鲸的数量'); ylabel('磷虾的数量'); 命令窗口输入：

第一问的三个M文件要分别在命令窗口输入：

17

一：euler(200,0.1) 二：impeuler(200,0.1) 三：rk4(200,0.1)

第二问的三个M文件要分别在命令窗口输入： 一：deuler(20,0.1) 二：dimpeuler(200,0.1) 三：drk4(20,0.1) （3） 运行结果: 第一问中的euler方法：

第一问中的改进的euler方法：

18

第一问中的经典的4级4阶Runge-Kutta法：

第二问中的euler方法：

19

第二问中的改进的euler方法：

20

第二问中的经典的4级4阶Runge-Kutta法：

（4） 结果分析：

三种算法的图像都会随着步长h的大小变化而变化。在第一问中，三个图像的差别不大；但在第二问中三种图像就有了明显的差别，而这种差别在于三种方法本身用于求解常微分问题的不同思路所致。

龚文胜 2003．7

21

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