数值分析实验报告

《数值分析》

实 验 指 导 书

徐州师范大学计算机科学与技术学院

2009年 11月

前 言

《数值分析实验指导书》根据数值分析实验教学的需要所编写，主要包括拉格朗日插值法、最小二乘法、数值积分、常微分方程初值问题数值解法、非线性方程求解、高斯消元法、数值计算方法应用等七个方面的内容。

前六个实验，都给出了实验目的、实验要求、实验原理和实验内容。实验原理主要介绍了算法的公式，并举例说明该方法的使用，通过这一部分，可以有效的将上课内容和实验内容对应联系起来。实验内容分为算法流程图、必做题和选做题三种形式。算法流程图在前两个实验给出，然后让学生逐步自己画流程图，锻炼了学生根据公式分析、书写流程图的能力，加深学生对算法的理解；学生通过编写必做题部分的程序、调试并运行，提高学生的编程能力；这两部分的内容是要求每一位学生都必须完成的。选做题是为基础较好的同学准备的，提高学生的进一步解决问题的能力。

本书不仅可以作为普通院校学生的数值分析上机实习指导教材，也可作为广大读者朋友掌握数值分析算法的实用教程。

由于时间仓促以及水平有限，错误与不足之处难免存在，敬请使用者给与批评指正！

目 录

实验一 拉格朗日插值法（2课时） .................................................................................................................. 1 一、实验目的 .................................................................................................................................................... 1 二、实验要求 .................................................................................................................................................... 1 三、实验原理 .................................................................................................................................................... 1 四、实验内容 .................................................................................................................................................... 2 实验二 最小二乘法（4课时） .......................................................................................................................... 5 一、实验目的 .................................................................................................................................................... 5 二、实验要求 .................................................................................................................................................... 5 三、实验原理 .................................................................................................................................................... 5 四、实验内容 .................................................................................................................................................... 6 实验三 数值积分（2课时） ............................................................................................................................ 10 一、实验目的 .................................................................................................................................................. 10 二、实验要求 .................................................................................................................................................. 10 三、实验原理 .................................................................................................................................................. 10 四、实验内容 ................................................................................................................................................... 11 实验四 常微分方程初值问题数值解法（4课时） ........................................................................................ 15 一、实验目的 .................................................................................................................................................. 15 二、实验要求 .................................................................................................................................................. 15 三、实验原理 .................................................................................................................................................. 15 四、实验内容 .................................................................................................................................................. 16 实验五 非线性方程求解（2课时） ................................................................................................................ 19 一、实验目的 .................................................................................................................................................. 19 二、实验要求 .................................................................................................................................................. 19 三、实验原理 .................................................................................................................................................. 19 四、实验内容 .................................................................................................................................................. 20 实验六 高斯消元法（2课时） ........................................................................................................................ 23 一、实验目的 .................................................................................................................................................. 23 二、实验要求 .................................................................................................................................................. 23 三、实验原理 .................................................................................................................................................. 23 四、实验内容 .................................................................................................................................................. 25 实验七 数值计算方法应用（6课时） ............................................................................................................ 27 一、实验目的 .................................................................................................................................................. 27 二、实验要求 .................................................................................................................................................. 27 三、实验内容 .................................................................................................................................................. 27 参考文献 ............................................................................................................................... 错误！未定义书签。

I

实验一 拉格朗日插值法（2课时）

一、实验目的

1.了解拉格朗日插值法的基本原理和方法。 2.掌握拉格朗日插值多项式的用法。 二、实验要求

1.用C语言编程实现拉格朗日插值。 2.进一步加深对拉格朗日插值法的理解。 三、实验原理

（一）拉格朗日插值公式

经过(n1)个点，(x0,y0),(x1,y1),(xn,yn)，构造一个n次多项式，形如：

pn(x)yklk(x)

k0n使得pn(xk)yk (k0,1,2,n)成立。

0(xx0)(xxk1)(xxk1)(xxn)其中lk(x)(xkx0)(xkxk1)(xkxk1)(xkxn)1xxkxxk

为插值基函数。 （二）例子

已知f(x)满足f(144)12,f(169)13,f(225)15作f(x)的二次拉格朗日插值多项式，并求f(175)的近似值。

解：设x0144,x1169,x2225，y012,y113,y215，则f(x)的二次拉格朗日插值基函数为：

0l(x)l(x)l(x)(xx1)(xx2)(x169)(x225) (x0x1)(x0x2)2025(xx0)(xx2)(x144)(x225)

(x1x0)(x1x2)1400(xx0)(xx1)(x144)(x169)(x2x0)(x2x1)4536 1 2因此f(x)的二次拉格朗日插值多项式为L2(x)y0l0(x)y1l1(x)y2l2(x)； 且f(175)L2(175) 12l0(175)13l1(175)15l2(175) 13.23015873。

1

四、实验内容

（一）算法流程图

如果程序的数据输入项（函数参数）为：插值节点及函数值，及待求点x的值；输出为待求点x对应的函数值，则程序流程图如图1-1。

输入xi,yi,及n,x,k=0 y = 0 输出k>n P =1 j=0,1…n k=j? P = P*(x-xj)/(xk-xj) y = y+P*yk k=k+1

图1-1 拉格朗日插值法算法流程图

（二）编程作业

编写拉格朗日插值法通用子程序，并用以下函数表来上机求f(0.15),f(0.31)。 x f (x) 0.0 0.39894 0.1 0.39695 0.195 0.39142 0.3 0.401 0.38138 0.36812 0.5 0.35206

代码如下：

/ * C语言源代码 */ #include main() { static float Lx[10],Ly[10];

2

int n,i,j; float x,y,p; printf(\"enter n=\"); scanf(\"%d\ // printf(\"enter xi\\n\"); for(i=1;i<=n;i++) scanf(\"%f\ printf(\"enter yi\\n\"); for(i=1;i<=n;i++) scanf(\"%f\ printf(\"enter x=\"); scanf(\"%f\ for( ) { ; for( ) { if(i!=j) ; // } ； } printf(\"y=%f\\n\} 请完成这个程序，并在两处注释处写上正确的注释。在执行程序时，如果求f(0.15)的值，那么屏幕上应出现如下内容： enter n=6 enter xi 0 0.1 0.195 0.3 0.401 0.5 enter yi 0.39894 0.39695 0.39142 0.38138 0.36812 0.35206

3

enter x=0.15 y=0.394473 （三）选做题

参考教材牛顿插值公式，编程实现用牛顿插值公式求上述条件下对应节点的函数值。

4

实验二 最小二乘法（4课时）

一、实验目的

1.了解最小二乘拟合的基本原理和方法，注意与插值方法的区别。 2.掌握最小二乘法。 二、实验要求

1.掌握用C语言作最小二乘多项式拟合的方法。 2.进一步加深对最小二乘法的理解。 三、实验原理

（一）最小二乘多项式拟合

已知数据对(xj,yj)(j1,2,,n)，求多项式

P(x)aixi(mn)

i02mi使得(a0,a1,,an)aixjyj为最小，这就是一个最小二乘问题。

j1i0nm

（二）法方程组为

mmxi i1...mxini1（三）最小二乘法计算步骤

xi1mi1mi.........2xim...n1xi...i1mxyia0i1i1mmn1axi 1  xiyii1i1.........manm2nnxxiiyii1i1nim  

用线性函数P(x)abx为例，拟合给定数据(xi,yi)(i1,2,,m)。 算法描述：

步骤1：输入m值，及(xi,yi)(i1,2,,m)。 步骤2：建立法方程组ATAXAY。 步骤3：解法方程组。 步骤4：输出P(x)abx。

5

四、实验内容

（一）算法流程图

1.算法整体流程图（如图2-1）

输入xi,yi(i=1,2,…,m)及m,n 生成中间矩阵C 生成法方程组的系数矩阵A=CTC 生成法方程组的右端向量b=CTY 解法方程中得ai(i=0,1,…,n) 输出ai(i=0,1,…,n)

图2-1 最小二乘法算法整体流程图

2.“生成中间矩阵C”算法流程图（如图2-2）

i=1,2,…,m C[i][1]←1 j=2,3,…,n+1 C[i][j]←x[i]*C[i][j-1]

图2-2 最小二乘法中“生成中间矩阵C”算法流程图

6

3.中间矩阵C的重要作用

1x11x2C1xm1x12则CTx1x1nx1x2xm222........................n21x2x2x2...mnx1y1nyx2,Y 2nxmym1xm2x2,nxmmmxTiACCi1mxini1myimi1xyiii1bCTYm2xiyii1mnxiyii1xixi1i1m2ixixi1i1mm23ixi1mn1ixi1mn2inxii1m...xin1,i1...m2n...xii1m...

（二）编程作业

测得铜导线在温度Ti(℃)时的电阻Ri()如下表，求电阻R与温度T的近似函数关系。

i Ti(℃) Ri() 0 19.1 76.30 1 25.0 77.80 2 30.1 79.25 3 36.0 80.80 4 40.0 82.35 5 45.1 6 50.0 83.90 85.10

【提示】在进行程序实现时，务必注意中间矩阵的作用，以及非奇次线性方程组求解问题！

为了实验的顺利完成，此处给出解非奇次线性方程组的高斯消元法的函数。请认真阅读并理解。

7

float gs(float a[20][20],float b[20],int n ) {int i,j,k,l; float s; k=1; while(k!=n+1) { if(a[k][k]!=0) { for(i=k+1;i<=n+1;i++) { a[i][k]=a[i][k]/a[k][k]; b[i]=b[i]-a[i][k]*b[k]; for(j=k+1;j<=n+1;j++) a[i][j]=a[i][j]-a[i][k]*a[k][j]; } } k=k+1; } for(k=n+1;k>=1;k--) { s=0; for(l=k+1;l<=n+1;l++) s=s+a[k][l]*b[l]; b[k]=(b[k]-s)/a[k][k]; } return 0; }

实验主程序如下（请加上必要的注释）。

8

9

实验三 数值积分（2课时）

一、实验目的

1.了解数值积分的基本原理和方法。 2.掌握复合梯形公式。

3.了解求积公式外推思想、Romberg公式及Romberg积分法。 二、实验要求

1.编写定步长复合梯形公式。 2.编写变步长复合梯形公式。

3.进一步加深对数值积分的理解。 三、实验原理

（一）定步长复合梯形公式 1.公式

将积分区间[a,b]n等分，分点为xiaih，i0,1,,n，其中h长。

ba称为积分步n2.例子

ban1hf(x)dxTnf(a)2f(xk)f(b)

2k14dx的近似值。用复合梯形公式求积分0（取8位小数，精确解为21x13.14159265）

11131T8{f(0)2[f()f()f()f()168482 537f()f()f()]f(1)}3.138988848（二）变步长复合梯形公式

1.公式

1banT2nTn22nk1bafa(2k1) 2n递推公式：

baf(a)f(b)T122k11babaTkTk1fa(2i1)2222ki12k(k1,2,3,)

2.例子

10

用递推公式求积分T146的近似值，使误差不超过。 10dx01x211f(0)f(1)142322111 T2T1f()3.12221113T4T2f()f()3.131176472444四、实验内容

（一）算法流程图

1.定步长复合梯形算法流程图

11

2.变步长复合梯形算法流程图 （二）编程作业

4求 0的近似值。

1x21（1）编写定步长复合梯形程序求解上式；

（2）编写变步长复合梯形程序求解上式，使误差不超过10-6。 【提示】请根据前面的算法流程图进行编写程序。

12

1.定步长复合梯形程序求解

13

2.变步长复合梯形程序求解

14

实验四 常微分方程初值问题数值解法（4课时）

一、实验目的

1.掌握Euler法和改进的Euler法公式的用法。 2.进一步加深对微分方程数值解的理解。

二、实验要求

1.编写欧拉法程序。

2.编写改进的欧拉法程序，学会用改进的欧拉公式来求解常微分方程初值问题。 三、实验原理

（一）常微分方程初值问题

(式1) 的数值解法，这也是科学与工程计算经常遇到的问题。 （二）欧拉法

求初值问题(式1)的一种最简单方法是将节点的导数代替，于是(式1)的方程可近似写成

用差商

(式2)

从出发(式2)右端，得到

称为解初值问题的Euler法。 （三）改进欧拉法

先用Euler法计算出

的近似

，将隐式梯形公式改为

的近似

，由(式2)求得

，一般写成

再将

代入

称为改进Euler法，它实际上是显式方法。即

（四）例子

2xyy，0x1y求初值问题的数值解，取步长h0.1。（精确解为

y(0)1

15

y(x)(12x)1/2）

yi11.1yi0.2xi/yi解：（1）利用欧拉法

y1,i0,1,2,90yi1yi0.05(K1K2)Ky2x/y1iii（2）利用改进欧拉法2(xi0.1)

Ky0.1Ki12yi0.1K1y01,i0,1,2,9计算结果如下： i xi Euler方法yi 改进Euler法yi 0 0 1 1 1 0.1 1.1 1.095909 2 0.2 1.191818 1.184096 3 0.3 1.277438 1.266201 4 0.4 1.358213 1.343360 5 0.5 1.435133 1.416402 6 0.6 1.508966 1.485956 7 0.7 1.580338 1.552515 8 0.8 1.649783 1.616476 9 0.9 1.717779 1.678168 10 1 1.784770 1.737869 精确解y(xi) 1 1.095445 1.183216 1.264991 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051 四、实验内容

（一）算法流程图

1.请根据欧拉公式，画出其算法流程图。

16

2.请根据改进欧拉公式，画出其算法流程图。

（二）编程作业

dyy,dx编写Euler法和改进的Euler法程序。求微分方程y(0)1x[0,1]在区间10等份

的近似解。

1.Euler法程序

17

2.改进的Euler法程序

在编写改进的Euler法程序时，有关输入输出部分，可参照以下屏幕上应出现内容。 please input a,b,h and a0 0 1 0.1 1

x=0.000000, y=1.000000 x=0.100000, y=0.905000 x=0.200000, y=0.819025 x=0.300000, y=0.741218 x=0.400000, y=0.670802 x=0.500000, y=0.607076 x=0.600000, y=0.549404 x=0.700000, y=0.497210 x=0.800000, y=0.449975 x=0.900000, y=0.407228 x=1.000000, y=0.368541 （三）选做题

使用梯形公式编写程序，解决上述例子问题。

18

实验五 非线性方程求解（2课时）

一、实验目的

1.了解求解非线性方程的解的常见方法。 2.编写牛顿迭代法程序求解非线性方程。 二、实验要求

1.设计牛顿迭代法算法，编写程序上机调试。 2.进一步加深对迭代法的理解。

三、实验原理

（一）牛顿迭代法

又称为牛顿-雷夫生方法（Newton-Raphson method），是一种在实数域和复数域上通过迭代计算求出非线性方程的数值解方法。方法的基本思路是利用一个根的猜测值x0做初始近似值，使用函数f(x)在x0处的泰勒级数展式的前两项做为函数f(x)的近似表达式。由于该表达式是一个线性函数，通过线性表达式替代方程f(x)= 0中的f(x)求得近似解 x1。即将方程f(x)= 0在x0处局部线性化计算出近似解x1，重复这一过程，将方程f(x)= 0在x1处局部线性化计算出x2，求得近似解

x2，„„。详细叙述如下：假设方程的解x*在x0 y 附近（x0是方程解x*的近似），函数f(x)在点x0 处的局部线化表达式为

f(x)f(x0)(xx0)f(x0)

由此得一次方程

f(x0)(xx0)f(x0)0

x 求解，得

f(x0)x1x0f(x0)

图5-1 牛顿迭代法示意图

如图5-1所示，x1比x0更接近于x*。该方法的几何意义是：用曲线上某点（x0，y0）的切线代替曲线，以该切线与x轴的交点（x1，0）作为曲线与x轴的交点（x*，0）的近似（所以牛顿迭代法又称为切线法）。设xn是方程解x*的近似，迭代格式

xn1xnf(xn) ( n = 0，1，2，„„) f(xn)就是著名的牛顿迭代公式，通过迭代计算实现逐次逼近方程的解。牛顿迭代法的最大优点是收敛速度快，具有二阶收敛。以著名的平方根算法为例，说明二阶收敛速度的意义。

（二）例子

已知21.4，求2等价于求方程f(x) = x2 – 2 = 0的解。由于f(x)2x。应用牛顿迭代法，得迭代计算格式

19

1（n = 0，1，2，„„） (xn2/xn)，

2取x0= 1.4为初值，迭代计算3次的数据列表如表5-1：

表5-1 牛顿迭代法数值实验 迭代次数 近似值 15位有效数 0 1.4 1.41421356237310 1 1.41428571428571 1.41421356237310 2 1.41421356421356 1.41421356237310 3 1.41421356237309 1.41421356237310 xn1误差 -1.42e-002 7.21e-005 1.84e-009 -2.22e-016 四、实验内容

（一）算法流程图

牛顿迭代法算法流程图

20

（二）编程作业

编写Newton迭代法通用子程序。实现方程f(x)=x6-x-1=0的满足精度要求的解。 要求求解过程中用一个变量I控制三种状态，其中： i=0表示求解满足给定精度的近似解； i=1表示f`(x0)=0，计算中断；

i=2表示迭代n次后精度要求仍不满足。 （说明：用教材P139例8）

21

（三）选做题

1.用简化牛顿法计算上述例题。 2.用牛顿下山法计算上述例题。

22

实验六 高斯消元法（2课时）

一、实验目的

1.掌握高斯选主元消去法公式的用法，适用范围及精确度。 2.通过高斯选主元消去法求矩阵方程的解，验证高斯消去法。

二、实验要求

1.写出高斯选主元消去法解线性方程组算法，编写程序上机调试出结果。 2.进一步加深对高斯消去法的理解。 三、实验原理

（一）高斯消元法

Gauss消去法就是将方程组

a11x1a12x2...a1nxnb1axax...axb2112222nn2 ................................................an1x1an2x2...annxnbn通过(n-1)步消元，将上述方程组转化为上三角方程组

再回代求此方程组的解。

（二）选列主元素高斯消去法

给定线性方程组Ax=b，记A(1)=A，b(1)=b，列主元Gauss消去法的具体过程如下： 首先在增广矩阵B(1)=(A(1),b(1))的第一列元素中，取

(1)(1)akmaxa为主元素，rkr1. 1i11in然后进行第一步消元得增广矩阵B(2)=(A(2),b(2))。 再在矩阵B(2)=(A(2),b(2))的第

二列元素中，取

(2)(2)ak为主元素，rkr2. 2maxai22in然后进行第二步消元得增广矩阵B(3)=(A(3),b(3))。按此方法继续进行下去, 经过n-1

步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n)).则方程组A(n)x=b(n)是与原方程组等价的上三角形方程组,可进行回代求解.

易证，只要|A|0,列主元Gauss消去法就可顺利进行。

（三）例子

采用4位十进制浮点计算，分别用顺序Gauss消去法和列主元Gauss消去法求解线

23

性方程组：

0.012x10.01x20.167x30.6781 x10.8334x25.91x312.13200x1200x4.2x981123方程组具有四位有效数字的精确解为 x1*=17.46，x2*=-45.76，x3*=5.546

解（1）用顺序Gauss消去法求解，消元过程为

0.01201.00032000.0120~000.0120~000.01000.833412000.1670.67815.91012.104.200981.08.01044.41 44531017981020.16700.67818.01044.41117510565171050.16700.67810.01000.100010314670.01000.10001030回代得： x3=5.546，x2=100.0，x1=-104.0

（2）用列主元Gauss消去法求解，消元过程为

0.01200.01000.16700.67811.0000.83345.91012.10320012004.200981.0选主元12004.200981.03200~1.0000.83345.91012.10 r1r30.01200.01000.16700.678112004.200981.03200~00.45845.90911.7900.551020.16700.67443200~003200~0012000.4584012000.45840981.05.90911.790.09610.53294.200

981.05.90911.790.09610.53294.200回代得：x3=5.545，x2=-45.77，x1=17.46

24

四、实验内容

（一）算法流程图

高斯列主元消去法N-S图

25

（二）编程作业

编写选列主元的高斯消去法。求出下列线性方程组Ax=b的解x。

321 16x1213x27

411x3

26

实验七 数值计算方法应用（6课时）

一、实验目的

初步具有使用数值计算方法解决实际问题的能力。 二、实验要求

在课后完成问题的分析，并编程上机求解。

三、实验内容

（一）使用课程有关“数值微分与数值积分”的知识上机解决以下问题： 1.贷款购房：某人从银行贷款购房，若他今年初贷款10 万元，月利率0.5%，每月还1000 元，建立差分方程计算他每年末欠银行多少钱，多少时间才能还清？如果要10 年还清，每月需还多少？

2.养老金基金：一老人60 岁时将养老金10 万元存入基金会，月利率0.4%，他每月取1000 元作生活费，建立差分方程计算他每岁末尚有多少钱？多少岁时将基金用完？如果60 岁时存入15万元，可以用到多少岁？

（二）使用课程有关“非线性方程及非线性方程组的解法”的知识上机解决以下问题： 1.某人以按揭方式贷款买了1 套价值20 万元的房子，首付了5 万元，每月还款1000元，15 年还清。问贷款利率是多少？

2.某人欲贷款50 万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15 年还清; 第二家银行开出的条件是每年还45000 元，20 年还清。从利率方面看，哪家银行较优惠（简单地假设年利率=月利率×12）？

实验一 拉格朗日插值法 #include void main() {

static float Lx[10],Ly[10]; int n,i,j;

float x,y=0,p;

printf(\"enter n=\");//输入点的个数 scanf(\"%d\

printf(\"enter xi\\n\");//输入点对应的x的值 for(i=1;i<=n;i++) scanf(\"%f\

printf(\"enter yi\\n\");//输入点对应的y的制 for(i=1;i<=n;i++) scanf(\"%f\

printf(\"enter x=\");//输入验证的x的值 scanf(\"%f\ for(i=1;i<=n;i++)

27

{

p=1;

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

if(i!=j)

p=p*(x-Lx[j])/(Lx[i]-Lx[j]); }

y=y+p*Ly[i]; }

printf(\"y=%f\\n\输入验证的x对应的y的值 }

实验二 最小二乘法 #include \"stdio.h\"

float gs(float a[20][20],float b[20],int n ) {

int i,j,k,l; float s; k=1;

while(k!=n+1) {

if(a[k][k]!=0) {

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

a[i][k]=a[i][k]/a[k][k]; b[i]=b[i]-a[i][k]*b[k]; for(j=k+1;j<=n+1;j++) a[i][j]=a[i][j]-a[i][k]*a[k][j]; } }

k=k+1; }

for(k=n+1;k>=1;k--) { s=0;

for(l=k+1;l<=n+1;l++) s=s+a[k][l]*b[l]; b[k]=(b[k]-s)/a[k][k]; }

return 0; }

28

int main() {

float a[20][20]={0.0};//定义a矩阵 float c[20][20];//定义c矩阵 float ct[20][20];//定义ct矩阵

float x[20];//定义数组用于存放x的数据 float y[20];//定义数组用于存放y的数据 float b[20]={0.0};//定义b矩阵

int i,j,k,m,n;

printf(\"输入所求函数的最高次数n:\\n\");//输入n（求线性的函数输入1。。） scanf(\"%d\

printf(\"输入测试数据的组数m:\\n\");//输入测试数据的组数 scanf(\"%d\

printf(\"输入x的测试数据%d个:\\n\输入x的测试数据m个 for(i=1;i<=m;i++) scanf(\"%f\

printf(\"输入y的测试数据%d个:\\n\输入y的测试数据m个 for(i=1;i<=m;i++) scanf(\"%f\

for(i=1;i<=m;i++)//c矩阵第一列赋值为1 c[i][1]=1.0;

//求C[][]

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

c[i][j]=x[i]*c[i][j-1];

//输出C[][]

printf(\"C矩阵如下：\\n\"); for(i=1;i<=m;i++) for(j=1;j<=n+1;j++) {

printf(\"%f \ if(j==n+1) printf(\"\\n\"); }

//求c的转置矩阵CT[][]

29

for(i=1;i<=m;i++) for(j=1;j<=n+1;j++) ct[j][i]=c[i][j];

//输出CT[][]

printf(\"CT矩阵如下：\\n\"); for(i=1;i<=n+1;i++) for(j=1;j<=m;j++) {

printf(\"%f \ if(j==m) printf(\"\\n\"); }

//求a[][]

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

a[i][j]+=ct[i][k]*c[k][j];

//输出a[][]

printf(\"a矩阵如下：\\n\"); for(i=1;i<=n+1;i++) for(j=1;j<=n+1;j++) {

printf(\"%f \ if(j==n+1) printf(\"\\n\"); }

//求b[]

for(i=1;i<=n+1;i++) for(k=1;k<=m;k++) b[i]+=ct[i][k]*y[k]; //输出b[]

printf(\"b矩阵如下：\\n\"); for(i=1;i<=n+1;i++) printf(\"%f \

gs(a,b,n);//调用高斯函数求方程组的解

30

printf(\"\\n\\n\");

printf(\"输出求得的函数的系数为:\\n\");

for(i=1;i<=n+1;i++)//输出求得的函数的系数 printf(\"a%d=%f \ printf(\"\\n\\n\");

return 0; }

实验三 数值积分 程序一

#include \"stdio.h\" float f(float x) {

return 4/(1+x*x); }

int main() {

int i,n;

float a,b,h,xi; float T=0.0;

printf(\"input a,b,n:\\n\");

scanf(\"%f%f%d\ T=T+f(a)+f(b); h=(b-a)/n; xi=a;

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

xi=xi+h; T=T+2*f(xi); }

T=T*h/2;

printf(\"求得的结果为：%f\\n\

return 0; }

程序二

#include \"stdio.h\" #include \"math.h\"

31

double f(double x) {

return 4.0/(1+x*x); }

int main() {

int n=1,k=1;

double a,b,T2n,Tn; printf(\"input a,b:\\n\"); scanf(\"%lf%lf\ Tn=(f(a)+f(b))*(b-a)/2.0;

T2n=Tn/2.0+(b-a)/2.0*f(a+(b-a)/2.0); printf(\"T1=%lf\\n\ printf(\"T2=%lf\\n\

while(fabs(T2n-Tn)>=0.000001) {

n=n*2; Tn=T2n; T2n=Tn/2.0;

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

T2n=T2n+f(a+(2*k-1)*(b-a)/(2*n))*(b-a)/(2*n); printf(\"T%d=%lf\\n\ }

printf(\"\\n满足条件的结果为：%lf\\n\

return 0; }

实验四 常微分数值求解问题 程序一

#include \"stdio.h\"

double f(double xi,double yi){

return -yi; }

int main() {

int n;

double a,b,y0;

printf(\"input a,b,n,y0:\\n\"); scanf(\"%lf\ scanf(\"%lf\

32

scanf(\"%d\ scanf(\"%lf\ double h=(b-a)/n; double xi=a; double yi=y0;

printf(\"xi\欧拉格式yi\\n\"); printf(\"%lf\%lf\\n\ int i;

double xi1,yi1; for(i=1;i<=n;i++) {

xi1=xi+h;

yi1=yi+h*f(xi,yi); xi=xi1; yi=yi1;

printf(\"%lf\%lf\\n\ }

return 0; }

程序二

#include \"stdio.h\"

double f(double xi,double yi){ return -yi; }

int main() { int n; double a,b,y0; printf(\"input a,b,n,y0:\\n\"); scanf(\"%lf\ scanf(\"%lf\ scanf(\"%d\ scanf(\"%lf\ double h=(b-a)/n; double xi=a; double yi=y0; printf(\"xi\改进的欧拉格式yi\\n\");

33

printf(\"%lf\%lf\\n\ int i; double xi1,yi1,yp,yc; for(i=1;i<=n;i++) { xi1=xi+h; yp=yi+h*f(xi,yi); yc=yi+h*f(xi1,yp); yi1=(1.0/2.0)*(yp+yc); xi=xi1; yi=yi1; printf(\"%lf\%lf\\n\ } return 0; }

实验5非线性方程求解 #include \"stdio.h\" #include \"math.h\" double f(double x) { return pow(x,6)-x-1; }

double f_(double x) { return 6*pow(x,5)-1; }

int main() { double x0,x,eps; int i,k,N; printf(\"请输入初始x0的值：\\n\"); scanf(\"%lf\ printf(\"迭代次数：\\n\"); scanf(\"%d\ printf(\"请输入误差e：\\n\"); scanf(\"%lf\

printf(\"k\xk\\n\"); printf(\"0\%.8lf\\n\

34

x=x0-f(x0)/f_(x0); if(!(f_(x0)==0)) { for(k=1;k<=N;k++) { printf(\"%d\%.8lf\\n\ if(fabs(x-x0)N) i=2; } else { i=1; } if(i==0) { printf(\"\\nx=%.8lf\\n\ printf(\"f(x)=%.30lf\\n\ printf(\"k=%d\\n\ } else if(i==1) { printf(\"给定的x0导致计算中断!!\\n\"); } else if(i==2) { printf(\"迭代N次后仍不满足精度要求!!\\n\"); }

return 0; }

实验六 高斯消元法

35

#include \"stdio.h\" #include \"math.h\" int main() { double A[10][10],b[10]; double m[10]; double eps=0.000001; int i,j,n,k,ik;

printf(\"输入矩阵的阶数n：\\n\"); scanf(\"%d\ printf(\"输入矩阵A:\\n\"); for(i=1;i<=n;i++) for(j=1;j<=n;j++) scanf(\"%lf\ printf(\"输入矩阵b:\\n\"); for(i=1;i<=n;i++) scanf(\"%lf\ double max; double temp; for(k=1;kif(k!=ik) { for(j=k;j<=n;j++) { temp=A[k][j]; A[k][j]=A[ik][j]; A[ik][j]=temp; } temp=b[k]; b[k]=b[ik]; b[ik]=temp;

36

} for(i=k+1;i<=n;i++) { m[i]=A[i][k]/A[k][k]; for(j=k;j<=n;j++) { A[i][j]=A[i][j]-A[k][j]*m[i]; } } for(j=k+1;j<=n;j++) { b[j]=b[j]-b[k]*m[j]; } } printf(\"A上三角矩阵：\\n\"); for(i=1;i<=n;i++) for(j=1;j<=n;j++) { printf(\"%12lf\ if(j==n) printf(\"\\n\"); } printf(\"b矩阵:\\n\"); for(i=1;i<=n;i++) { printf(\"%12lf\} b[n]=b[n]/A[n][n]; for(i=n-1;i>=1;i--) { for(j=n;j>i;j--) { b[i]=(b[i]-A[i][j]*b[j]); } b[i]=b[i]/A[i][i]; }

37

printf(\"\\n求得得结果为：\\n\"); for(i=1;i<=n;i++) {

printf(\"x%d=%lf\\ } return 0;

}

38