计算方法 用欧拉预估-校正法求初值问题

《计算方法》实验指导书

《计算方法》实验指导书

实验1 方程求根

一、实验目的

1． 通过对二分法、牛顿法、割线法作编程练习，进一步体会它们各自不同的特点； 2． 了解二分法，切线法，割线法。

3． 能熟练运用二分法，牛顿法进行方程求根

4． 通过上机调试运行，对方程求根的几种方法程序进行改进。

二、实验要求

1． 上机前作好充分准备，包括复习编程所需要的语言工具。 2． 上机时要遵守实验室的规章制度，爱护实验设备。

3． 记录调试过程及结果，记录并比较与手工运算结果的异同。 4． 程序调试完后，须由实验辅导教师在机器上检查运行结果。 5． 给出本章实验单元的实验报告。

三、实验环境、设备

1． 硬件设备：IBM PC以上计算

机，有硬盘和一个软驱、单机和网络环境均可。

2． 软件环境： C语言运行环境。

四、实验原理、方法 二分算法计算步骤：

（1）输入有根区间的端点a、b及预先给定的精度ε；

（2）计算中点x=(a+b)/2；

（3）若f(x)f(b)<0，则a=x，转向下一步；否则b=x，转向下一步； （4）若b-a<ε，则输出方程满足精度要求的根x，结束；否则转向步骤（2）。 迭代法：

256

《计算方法》实验指导书

开始 输入 x0, ε,N k=1 开始 输入x0, ?,N k=1 f'(x0)=0? 是 k=k+1 x1=?(x0) 是 否 x?x?10f(x)0f'(x)0x0=x1 x1?x0≤ε k=k+1 x0=x1 是 输出奇异标志 是 否 是 k

输出近似根x1

图 2.3 迭代法框图

牛顿法：

牛顿迭代法是一种逐步线性化方法，即将非线性方程f(x)=0的求根问题归结为计算一系列线性方程的根。

设xk是方程f(x)=0的一个近似根，将f(x)在xk处作一阶泰勒展开，即

f(x)≈f(xk)+f′(xk)(x- xk)

于是得到如下的近似方程

f(xk)+f′(xk)(x- xk)=0 （2.7）

设f′(xk)≠0，则式（2.7）的解为

x?xk?取x作为原方程的新的近似根xk+1，即令

f(xk) 'f(xk)xk?1?xk?f(xk) k=0,1,2, … （2.8）

f'(xk)则称式（2.8）为牛顿迭代公式。用牛顿迭代公式（2.8）求方程近似根的方法称为牛顿迭代法，简称牛顿法，又称切线法。

五、实验内容

1． 以方程：x-0.2x-0.2x-1.2=0为例，编写程序求方程的根

257

3

2

《计算方法》实验指导书

2． 编写二分法、迭代法、牛顿法程序，分析运行结果。 3． 对用这两种方法求解出的根进行对比分析

六、实验步骤

1． 根据实验题目，给出题目的C程序。 2． 上机输入和调试自己所编的程序。 3． 上机结束后，应整理出实验报告。

七、实验报告要求及记录、格式

按金陵科技学院《实验报告（工科）》格式填写

附1：牛顿法程序核心部分：

for(i=0;i

{printf(\

x1=x0-f(x0)/f1(x0); /*牛顿迭代*/

if(fabs(x1-x0)

{printf(\满足精度，输出近似根*/ return; } x0=x1; }

258

《计算方法》实验指导书

实验2 线性方程组数值解法

一、实验目的

1．掌握方程组的解法，迭代法及其收敛性。

2．能熟练掌握高斯消去法，列主元高斯消去法，三角分解法。 3．掌握雅可比迭代法，高斯=赛德尔迭代求线性方程组的解。

二、实验要求

1．上机前作好充分准备，比较不用的方法解决相同问题的不同。 2．上机时要遵守实验室的规章制度，爱护实验设备。 3．记录调试过程及结果，记录并比较与手工

开始 运算结果的异同。

4．程序调试完后，须由实验辅导教师在机器

输入aij,bi和上检查运行结果。

方程阶数n 5．给出本章实验单元的实验报告。

三、实验设备、环境

1．硬件设备：IBM PC以上计算机，有硬盘和一个软驱、单机和网络环境均可。 2．软件环境： C语言运行环境。

是 K=1 列主元 四、实验原理、方法

1、高斯消去法：

1）计算步骤

(1)输入方程组的阶数n，系数矩阵Ａ和右端常数矩阵b

(2)消元过程：设a计算

(k)k?mik?aik/akk?(k?1)(k)(k)k+2，…，?aij?aij?mikakj i,j=k+1，

?(k?1)(k)(k)b?b?mbiiikk?(k)kk(k)akk?0? 否 打印奇标志 aik=aik/akk i=k+1,k+2,…,n 结束 aij= aij -aik*akj bi =bi -aik*bk i,j=k+1,k+2,…,n 否 k=n-1? 是 k=k+1 ?0，对k=1，2，…，n-1

bn=bn/ann 1bi?(bi?anj?i?1n

(3)回代过程

(n)?bn?xn?a(n)?nn??x?(b(i)?ii???ab) ijjni=n-1,n-2,…,1 输出 (i)(i)xj)/aii,i?n?1,?,2,1?aijnb1,b2,…,bn 结束 高斯消去法框图

j?i?1

(４)输出方程组的解

259

《计算方法》实验指导书

2、列主元高斯消去法

(1)、输入方程组的阶数n，系数矩阵Ａ和右端常数矩阵b (2)、列主元素：对k=1，２，…，n-1，选出akk?(k?1)(k?1)(k?1)(k?1)中绝对值最大的元素am,k，,ak?1,k,?,ank?k行和m行交换后，再作第k步消元操作。

(3)、消元过程：对k=1，2，…，n-1计算

?(k)(?m?ak)ikik/akk?a(k?1)(k)m(k)?ij?aij?ikakj ?b(k?1)i?b(k)i?mikb(k)k（i,j=k+1，k+2，…，n） (4)、回代过程

?(??xbn)nn?(?an)nn

?)?x?(b(i)?x(i)?iiijj)/aii(i?n?1,?,2,1)j?na(i?i?1(5)、输出方程组的解

3、三角分解法：

（1）根据方程组得到增广矩阵 （2）对j=1，2，…，n计算u1j?a1j

对i=2，3，…，n计算lai1i1?u

11（3）对k=1，2，…，n

k?1 a.对j=k,k+1,…,n+1计算ukj?akj??lkquqjq?1

b.对i=k+1，k+2，…，n计算lik?(aik??k?1liquqk)/ukkq?1

（4）回代计算解

xn?yn/unn，对k=n-1，n-2，…，2，1计算

nxk?(yk?kqxq)/ukk

q?u?k?1

260

从主程序来 d=akk L=k i=k+1 否 aik?d是 d=aik L=i 否 i=n? i=i+1 是 是 d=0? 是 否 打印奇标志 L=k? 否结束 t=aij,aij=akj,akj=t j=k,k+1,…,n t=bk,bk=bL,bL=t 返回主程序 列主元框图

对

《计算方法》实验指导书

???f[x0,x1,?,xn](x?x0)(x?x1)?(x?xn?1)

（4）、输出函数值。

五、实验内容

1．已知函数表 x 1.615 1.634 1.702 1.828 1.921 Y=f(x) 2.41450 2.46459 2.65271 3.03035 3.34066 （1）用拉格朗日插值（二次、四次）计算f(1.682)和f(1.813)的近似值。

（2）构造出均差表，并利用牛顿（均差）插值多项式计算f(1.682)和f(1.813)的近似值。

（3）分析并比较两种算法得到的近似值的精度。

2．编写拉格朗日和牛顿插值算法程序，分析运行结果。

六、实验步骤

1． 根据实验题目，给出解决问题的程序代码。 2． 上机输入和调试自己所编的程序。 3． 上机结束后，应整理出实验报告。

七、实验报告要求及记录、格式

按金陵科技学院《实验报告（工科）》格式填写

附1：拉格朗日插值法核心程序：

for(i=0;i

for(j=0;j

if(j!=i) f[i]*=(xx-x[j])/(x[i]-x[j]); yy+=f[i]; }

附2：牛顿插值法核心程序：

for(I=1;I<=N;I++) /*计算牛顿插值函数的值*/ { f[0]=y[I];

for(j=0;j

b=y[N];

for(I=N-1;I>=0;I--) b=b*(xx-x[I])+y[I]; /*计算函数值*/

266

开始 输入(xi,yi),N i=0,1,2,…,N k=1 f(i+l)=f(i)?u(i)x(k)?x(i) i=0,1,…,k-1 否 k=N? y(k)=f(k) k=k+1 是 输出x1,x2 N1=f(N);N2=f(N) N1=N1*(x1-x(i))+y(i) N2=N2*(x2-x(i))+y(i) i=0,1,…,k-1 输出N1;N2 结束 牛顿插值法框图

《计算方法》实验指导书

实验4 曲线拟合

一、实验目的

1．掌握最小二乘原理，正规方程组，超定方程组概念。 2．掌握用最小二乘法拟合曲线，超定方程级的最小二乘解。 3．掌握用最小二乘法拟合曲线。

二、实验要求

1．上机前作好充分准备，复习最小二乘拟合方法。 2．上机时要遵守实验室的规章制度，爱护实验设备。

3．记录调试过程及结果，记录并比较与手工运算结果的异同。 4．程序调试完后，须由实验辅导教师在机器上检查运行结果。 5．给出本章实验单元的实验报告。

三、实验设备、环境

1．硬件设备：IBM PC以上计算机，有硬盘和一个软驱、单机和网络环境均可。 2．软件环境： C语言运行环境。

四、实验原理、方法

最小二乘法：

对于给定的线性方程组

Ax=b 式中

?b1??x1??a11 a12 ? a1n??b??x??a a? a?221222n?，b=??，x=?2? A=?????????????????????ba a? a?m??xn?mn??m1m2当m>n时，称为矛盾方程组，又称超定方程组。对于这种方程组有m个方程，而只有比

m小的n个变量，即方程的个数超过未知量的个数，这种方程组一般来说是没有解的。我们转而寻求在某种意义下的近似解。如果这组近似解对于矛盾方程组中的每个方程式的误差的平方和为最小，即

?(x1,x2,?,xn)?????(?aijxj?bi)2

2ii?1i?1j?1mmn为最小值时，就可以认为该组近似值为矛盾方程组的近似解。 这种近似解不是指对精确解的近似（因为精确解并不存在），而是指寻求各未知数的一组取值，使方程式（5.1）中各方程式近似相等，这就是最小二乘法的基本思想。

五、实验内容

1．设有如下实验数据

267

《计算方法》实验指导书

x y

1.36 1.49 1.73 1.81 1.95 2.16 2.28 2.48

14.094 15.069 16.844 17.378 18.435 19.949 20.963 22.495

试用最小二乘法分别求一次及二次多项式曲线拟合以上数据。 2．编写程序，分析运行结果。

六、实验步骤

1． 根据实验题目，给出解决问题的程序代码。 2． 上机输入和调试自己所编的程序。 3． 上机结束后，应整理出实验报告。

七、实验报告要求及记录、格式

按金陵科技学院《实验报告（工科）》格式填写

附：抛物函数拟合源程序

/********************************************************/ /* 最小二乘法-拟合抛物函数 */ /********************************************************/ #include #include

#define Max_N 25 /*最大数据点的个数*/ #define M 3 /*正规方程组的阶数*/ /*列主元高斯消去法求解线性方程组*/ void ColPivot(float A[M][M],float B[],int n) { int i,j,k,m_i; float m_x,temp; for(i=0;i

j=i+1; m_i=i; m_x=fabs(A[i][i]); for(;j

if(fabs(A[j][i]>m_x)) /*找主元素*/ {m_i=j;

m_x=fabs(A[j][i]); }

if(i

{ temp=A[i][j]; A[i][j]=A[m_i][j]; A[m_i][j]=temp; } }

/*消元*/

for(j=i+1;j

268

《计算方法》实验指导书

B[j]+=B[i]*temp; for(k=i;k

A[j][k]+=A[i][k]*temp; } } }

main() {int i,j,k,n;

float x[Max_N],y[Max_N],b[M],a[M][M],c[M]; printf(\输入点数n*/ do

{ scanf(\ if(n>Max_N)

printf(\ }

while(n>Max_N||n<=0);

printf(\ for(i=0;i

a[i][j]=a[i][j]+pow(x[k],i+j); b[i]=b[i]+pow(x[k],i)*y[k]; } } }

/* 输出正规方程组 for(i=0;i

ColPivot(a,b,M); /*调用列主元消去法函数计算方程组的解*/

c[M-1]=b[M-1]/a[M-1][M-1]; /*解方程*/ for(i=M-2;i>=0;i--) { c[i]=b[i];

for(j=i+1;j

269

《计算方法》实验指导书

c[i]-=a[i][j]*c[j]; c[i]/=a[i][i]; }

printf(\ /*输出方程组的解*/ for(i=0;i

printf(\

printf(\：y=%f+(%f)x+(%f)x2\ getch(); }

/*---------------------------------------- End of file ------------------------------------*/ /*程序输入输出

Please input n value:16 Input x[i],i=0,...N-1:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Input y[i],i=0,...N-1:

4 6.4 8 8.8 9.22 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.60 Solve is :

c[0]=4.387500 c[1]=1.065962 c[2]=-0.044466

Result:y=4.387500+1.065962x+(-0.044466)x2

*/

270

《计算方法》实验指导书

实验5 数值积分与数值微分

一、实验目的

1．掌握插值型求积公式，梯形公式、辛卜生公式、柯特斯公式。 2．掌握高斯求积公式，数值微分的中点公式。 3．掌握龙贝格求积公式。

二、实验要求

1．上机前作好充分准备，理解求积公式，梯形公式的概念。

2．上机时要遵守实验室的规章制度，爱护实验设备。

3．记录调试过程及结果，记录并比较与手工运算结果的异同。

4．程序调试完后，须由实验辅导教师在机器上检查运行结果。

5．给出本章实验单元的实验报告。

三、实验设备、环境

1．硬件设备：IBM PC以上计算机，有硬盘和一个软驱、单机和网络环境均可。 2．软件环境： C语言运行环境。

四、实验原理、方法

a) 梯形求积公式：

复化梯形公式计算步骤

(1)、令h=(b-a)/N，T=0（h为等分数，T为存

放积分值的变量）

(2)、对k=1,2,…,N计算

T=T+f(a+kh)

(3)、T=h[f(a)+2T+f(b)]/2 b) 辛卜生求积公式

复化辛卜生公式计算步骤

(1)、令h=(b-a)/N，s1=f(a+h/2)，s2=0 (2)、对k=1,2,…,N-1计算

s1=s1+f(a+kh+h/2)，s2=s2+f(a+kh) (3)、s=h[f(a)+4s1+2s2+f(b)]/2。

c) 龙贝格求积公式

算法流程见：龙贝格求积公式流程图

271

开始 定义f(x) 输入a,b,N h=(b-a)/N T=0 对于(k=1,2,…,N-1) T=T+f(a+kh) T=h[f(a)+2T+f(b)]/2 结束 复化梯形公式流程图 开始 定义f(x) 输入a,b,N h=(b-a)/N,x=a+h/2, s1=f(x),s2=0 s1=s1+f(a+kh+h/2) s2=s2+f(a+kh) (k=1,2,…,N-1) s2=h[f(a)+4s1+2s2+f(b)]/6 结束 复化辛卜生公式流程图

《计算方法》实验指导书

五、实验内容

1．分别编写出梯形公式、辛卜生求积公式、龙贝格求积公式算法程序。 2．分别用上述算法程序求I??1dx01?x的积分值，并分析每种算法程序求得的结果为什么不同。

六、实验步骤

1． 根据实验题目，给出解决问题的程序代码。2． 上机输入和调试自己所编的程序。 3． 上机结束后，应整理出实验报告。

七、实验报告要求及记录、格式

按金陵科技学院《实验报告（工科）》格式填写

附：龙贝格求积分法源程序

/**************************************/ /* 用龙贝格求积公式求定积分的解 */ /**************************************/ #include #include

#define epsilon 0.0001 /*精度要求，精度不能太高，否则不能得出结果*/ /*求积函数f(x)*/ float f(float x) {return(sin(x)); } /*梯形公式*/

float Romberg(float a,float b) {int k=1;

float S,x,T1,T2,S1,S2,C1,C2,R1,R2,h=b-a; T1=h*(f(a)+f(b))/2; while(1) {S=0; x=a+h/2; do { S+=f(x); x+=h; }while(x

if(fabs(T2-T1)

272

开始 输入a，b，ε h=b-a，k=1 T1=h[f(a)+f(b)]/2 S=0，x=a+h/2 S=S+f(x)，x=x+h 是 x

《计算方法》实验指导书

if(k==1) { T1=T2;S1=S2;h/=2;k+=1; continue;} C2=S2+(S2-S1)/15.0;

if(k==2) { C1=C2;T1=T2;S1=S2;h/=2;k+=1; continue;} R2=C2+(C2-C1)/63.0;

if(k==3) { R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1; continue;} if(fabs(S2-S1)

printf(“\\nInput the begin and end:”); /*输入积分区间*/ scanf(“%f%f”,&a,&b);

S=Romberg(a,b); /*调用龙贝格算法函数*/ printf(“Solve is: %f”,S); getch(); }

/*--------------------------- End of file ----------------------------*/ /*

对函数f(x)=sin(x)在积分区间[1,2]内求定积分。 程序输入输出

Input the begin and end: 1.0 2.0 Solve is: 0.956449 */

273

《计算方法》实验指导书

实习6 常微分方程数值解法

一、实验目的

1、通过本实验，充分理解常微分方程的初值问题的有关方法和理论理论；

2、通过实际计算体会各种解法的功能、优缺点及适用场合，会选取适当的求解方法。

二、实验要求

1、了解各种解决常微分方程初值问题的方法，并比较各种方法的不同之处，充分理解各种方法的特点及用途；

2、针对后面的练习题目进行上机计算；

3、在充分理解各种方法的基础上，会编写各种方法的程序； 4、考虑其他龙格-库塔公式的程序编写。

三、实验设备、环境

1．硬件设备：IBM PC以上计算机，有硬盘和一个软驱、单机和网络环境均可。 2．软件环境： C语言运行环境。

开始 输入x0,y0,h,N 四、实验原理、方法

1、 改进欧拉公式算法框图见右侧流程图 2、 龙格-库塔公式（标准四阶龙格-库塔公式）

开始 输入x0,y0,h,N n=n+1 x0=x1 y0=y1 n=1 x1= x0 +h yp= y0 +hf(x0,y0) yc= y0 +hf(x1,yp) y1=(yp+yc)/2 n=1 输出x1,y1 x1=x0+h k1=f(x0,y0),k2=f(x0+h/2,y0+hk1/2) k3=f(x0+h/2,y0+hk2/2),k4=f(x1,y0+hk3) y1=y0+h(k1+2k2+2k3+k4)/6 否 n=N? 是 结束 输出x1,y1 否

改进欧拉法流程图

n=n+1 x0=x1,y0=y1

n=N? 是 结束 龙格-库塔法流程图

274

《计算方法》实验指导书

五、实验内容

1、用欧拉方式和改进欧拉法求初值问题 y'=x*y; y(3.0)=1; 3.0<=x<=3.6 取h=0.2(即n=3) 2、龙格-库塔公式求初值问题 y'=x/y; y(2.0)=1; 2.0<=x<=2.6 取h=0.2(即n=3)。 3、编写欧拉方式和改进欧拉法算法程序，调试四阶龙格-库塔公式程序。

六、实验步骤

4． 根据实验题目，给出解决问题的程序代码。 5． 上机输入和调试自己所编的程序。 6． 上机结束后，应整理出实验报告。

七、实验报告要求及记录、格式

按金陵科技学院《实验报告（工科）》格式填写

附：龙格-库塔公式（标准四阶龙格-库塔公式）程序

例：求初值问题 y'=x/y; y(2.0)=1; 2.0<=x<=2.6 取h=0.2(即n=3)。 /****************************************************/ /* 标准四阶龙格-库塔公式 */ /****************************************************/ #include #include /*y1=f(x,y)*/ float f(float x,float y) {return(x/y);}

/*四阶龙格-库塔公式函数*/

void Runge_Kutta(float x0,float xn,float y0,int n) {int i;

float x=x0,y=y0,k1,k2,k3,k4,h=(xn-x0)/n; printf(\ for(i=1;i<=n;i++) {k1=f(x,y);

k2=f(x+h/2,y+h*k1/2); k3=f(x+h/2,y+h*k2/2); k4=f(x+h,y+h*k3);

y=y+h*(k1+2*k2+2*k3+k4)/6; x=x0+i*h;

printf(\ } } main() { int i,n; float h,x0,xn,y0;

printf(\ /*输入x的区间*/ scanf(\

printf(\ /*输入x0处的y的值y0*/

275

《计算方法》实验指导书

scanf(\ do

{ printf(\ /*输入区间的等分数*/ scanf(\ }

while(n<=1);

Runge_Kutta(x0,xn,y0,n); getch(); }

/*----------------------- End of file -------------------------------*/ /* 程序输入输出

Input the begin and end of x: 2.0 2.6 Input the y value at 2.000000: 1

Input n value[divide (2.000000,2.600000)]: 3 x[0]=2.000000 y[0]=1.000000 x[1]=2.200000 y[1]=1.356505 x[2]=2.400000 y[2]=1.661361 x[3]=2.600000 y[3]=1.939104 */

276

《计算方法》实验指导书

scanf(\ do

{ printf(\ /*输入区间的等分数*/ scanf(\ }

while(n<=1);

Runge_Kutta(x0,xn,y0,n); getch(); }

/*----------------------- End of file -------------------------------*/ /* 程序输入输出

Input the begin and end of x: 2.0 2.6 Input the y value at 2.000000: 1

Input n value[divide (2.000000,2.600000)]: 3 x[0]=2.000000 y[0]=1.000000 x[1]=2.200000 y[1]=1.356505 x[2]=2.400000 y[2]=1.661361 x[3]=2.600000 y[3]=1.939104 */

276

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