本科生课程设计报告实习课程数值分析
学院名称管理科学学院
专业名称信息与计算科学
学生姓名
学生学号
指导教师
实验地点
实验成绩
二〇一六年六月二〇一六年六
摘要
常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.
关键词：数值解法；常微分
1. 实验内容和要求
常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。
要求：分别取步长h = ，，，采用改进的Euler 方法、4阶经典龙格-库塔R －K 方法和4阶Adams 预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。其中多步法需要的初值由经典R-K 法提供。运算时，取足以表示计算精度的有效位。
2. 算法说明
函数及变量说明
表1 函数及变量定义
1、欧拉方法：
1()()(,())i i i i y x y x hf x y x +=+ （i=0，1，2，3，......n-1）
(0)y a
= (其中a 为初值)
2、改进欧拉方法：
~
1~111()()(,())
()()[(,())(,())]
2
(0)i i i i i i i i i i y x y x hf x y x h
y x y x f x y x f x y x y a ++++=+=++=（i=0，1，2......n-1） (其中a 为初值)
3、经典K-R 方法： 1
1213
243()6
(,)(,)22(,)22
(,)
i i i i i i i i i i h y y K f x y h h
K f x y K h h K f x y K K f x h y hK +?
=+??
=???=++??
?
=++??=++??
4、4阶adams 预测-校正方法 预测： 校正：
Adsms 内插外插公式联合使用称为Adams 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。计算时需要注意以下两点： 1、外插公式为显式，内插公式为隐式。故用内插外插公式时需要进行迭代。
2、这种预测-校正法是四步法，计算Yn+1时，不但用到前一步信息，而且要用到更前三步信息1-n f ，2-n f 3-n f ，因此它不是自动开始的，实际计算时必须借助某种单步法，用Runge-Kutta 格式为它提供初始值
3
21,,y y y ，依据局部截断误差公式得：
进一步将Adams 预测-校正系统加工成下列方案：
运用上述方案计算1n y +时，要先一步的信息n n n
c p y -,,y '
n 和更前一步的信息1-n y
因此在计算机之前必须给出初值1y 和11c -p ，1y 可用其他单步法来计算，11c -p 则一般令他为零。
3. 源程序
#include<>
#include<>
#include<>
#include<>
//微分方程
double f(double x,double y)
{
return (-y+cos(2*x)-2*sin(2*x)+2*x*exp(-x));
//return x*pow(y,-2)*;
}
//原函数
double p_f(double x)
{
return x*x*exp(-x)+cos(2*x);
//return pow+pow(x,2),;
}
//求精确解
double* Exact(double a,double b,double h)
{
int i;
double c = (b-a)/h;
double *y = new double [c+1];
for(i=0;i<=c;i++)
y[i]= p_f(a+i*h);
return y;
}
//欧拉方法
double* Euler(double a,double b,double h,double y0) {
int i;
double c = (b-a)/h;
double *y = new double [c+1];
y[0]=y0;
for(i=1;i<=c;i++)
{
y[i]= y[i-1]+ h* f(a+(i-1)*h,y[i-1]);
//printf("%f\n",f(a+(i-1)*h,y[i-1]));
}
return y;
}
//改进欧拉方法
double* Euler_Pro(double a,double b,double h,double y0) {
int i;
double c = (b-a)/h ,yb;
double *y = new double [c+1];
y[0]=y0;
for(i=1;i<=c;i++)
{
yb=y[i-1] +h* f(a+(i-1)*h,y[i-1]);
y[i] = y[i-1] + *h*( f(a+(i-1)*h,y[i-1]) + f(a+i*h,yb) );
}
return y;
}
//经典K-R方法
double* K_R(double a,double b,double h,double y0)
{
double K1,K2,K3,K4,x;
int i;
double c = (b-a)/h ;
double *y = new double [c+1];
y[0]=y0;
for(i=1;i<=c;i++)
{
x=a+(i-1)*h;
K1=f(x,y[i-1]);
K2=f(x+*h,y[i-1]+*h*K1);
K3=f(x+*h,y[i-1]+*h*K2);
K4=f(x+h,y[i-1]+h*K3);
y[i]=y[i-1] + h*( K1 + 2*K2 + 2*K3 + K4)/6;
}
return y;
}
//4阶Adams预测-校正方法
double **Adams(double a,double b,double h,double y0)
{
int i;
double *y;
double count = (b-a)/h;
double dy[100]={0},x[4]={0};
double p_0=0,p_1=0,c=0;
double **r = new double *[count+1];//动态初始化，储存预测值和校值
for(i=0;i<=count;i++)
r[i] = new double [2];
if(r==NULL)
{
printf("内存分配失败\n");
exit(0);
}
for(i=0;i { r[i][0]=0; r[i][1]=0; } y=K_R(a,b,h,y0); for(i=0;i<4;i++) { x[i]=a+h*i; dy[i]=f(x[i],y[i]); r[i][0]=y[i]; } i=3; while(x[3] { x[3] = x[3]+h; p_1=y[3]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24; //预估 c=p_1+251*(c-p_0)/270; //修正 r[i][0]=c; c=f(x[3],c); //求f c=y[3]+h*(9*c+19*dy[3]-5*dy[2]+dy[1])/24; //校正 y[3]=c-19*(c-p_1)/270; //修正 r[i][1] = y[3]; dy[0]=dy[1]; dy[1]=dy[2]; dy[2]=dy[3]; dy[3]=f(x[3],y[3]); p_0=p_1; i++; } return r; } void main() { int i , selet; double a=0, b=2, h=, y0=1; while(1) { printf("请输入下限，上限，步长，初值：(用空格隔开)\n"); scanf("%lf %lf %lf %lf",&a,&b,&h,&y0); double c = (b-a)/h; double *yx,*y; yx=Exact(a,b,h); while(1) { printf("1、欧拉 2、改进欧拉 3、经典K-R4、4阶Adams 预测-校正 5、修改数据 6、退出\n"); scanf("%d",&selet); system("cls"); if(selet==1) { y=Euler(a,b,h,y0); printf("欧拉方法：\n"); } if(selet==2) { y=Euler_Pro(a,b,h,y0); printf("改进欧拉方法：\n"); } if(selet==3) { y=K_R(a,b,h,y0); printf("经典K-R方法(4阶龙格-库塔方法)：\n"); } if(selet==4) { double **r; r=Adams(a,b,h,y0); printf("4阶Adams预测-校正方法:\n"); printf("x值预估值校正值误差\n"); printf("% %lf - -\n",0,r[0][0]); printf("% %lf - -\n",,r[1][0]); printf("% %lf - -\n",,r[2][0]); for(i=3;i<=10;i++) printf("% %f %f %\n",a+i*h,r[i][0],r[i][1],fabs( r[i][0]-r[i][1])); } if(selet==1||selet==2||selet==3) { printf("x值计算值准确值误差\n"); for(i=0;i<=10;i++) printf("% %lf %lf %\n",(a+i)*h,y[i],yx[i],fabs(yx[ i]-y[i])); } if(selet==5) break; if(selet==6) exit(1); continue; } } } 4.实验结果 1、改进欧拉方法 2、经典K-R方法 3、4阶adams预测-校正法 5.对比分析 通过表格可以看出： 1，在步长相同的时候，结果准确度经典K-R>改进欧拉>欧拉方法>4阶adams预测校正方法 2，同一个方法的情况下，步长越小结果准确度越高； 参考文献： [1] 朱建新，李有法，数值计算方法(第三版)，高等教育出版社； [2]Adam预测-校正系统 《计算机数学基础》数值部分第五单元辅导 14 常微分方程的数值解法 一.重点内容 1. 欧拉公式： ）心知1）a 儿+1 =儿 + hfg ，儿） m 1、 伙=0丄2,…川一 1） I 无=x Q +kh 局部截断误差是0（*）。 2. 改进欧拉公式： 预报一校正公式： 预报值 _v*+1 =儿+ hf （x k ,儿） - h - 校正值 y M = y k +-［f （x kt y k ） + /（x A+1, y M ）］ 即 儿+1 =儿+ ￡ "（忑'儿）+心+「儿+ hfg ,儿））］ 或表成平均的形式： 儿=儿+ hfg ,儿） '儿=儿+"（无+】，儿） +K ） 改进欧拉法的局部截断误差是0（2） 3. 龙格一库塔法 二阶龙格一库塔法的局部截断误差是0（爪） 三阶龙格一库塔法的局部截断误差是0（护） 四阶龙格F 塔法公式：儿计=儿+ 2（匕+ 2心+ 2? + ?） 四阶龙格一库塔法的局部截断误差是0（爪）。 二实例 y' = — y — xv f2（0 < x < 0.6） 例1用欧拉法解初值问题｛ ' ? -取步长/匸02计算过程保留 b （o ）= 1 4位小数。 解/i=0.2. f （x ）= —y —xy 2<,首先建立欧拉迭代格式 y*+i =儿+ hf g,y k ） = y k -hy k -hx k y ； =0.2 儿（4 - x k y k ）（k = 0,1,2） K 2=f(x n +^h, yk+-hK\)t gg+舟人，>'n +y/?A3）； 当k=0, xi=0.2 时，已知x（）=0,y（）=l,有 y(0?2)今i=0?2X l(4-0X 1)=0.8000 当k=\. M=0?4时，已知“=0?2」尸0?8,有 y(0?4)今2=0.2 X 0.8X(4-0.2X0.8)=0.614 4 当k=2, xs=0.6 时，已知x2=0.4,y2=0.6144,有 y(0?6)今3=0.2 X0.6144X (4-0.4 X 0.4613)=0.8000 「J, ,2 ?_ ZX 例2用欧拉预报一校正公式求解初值问题\y + v +V sinx=,取步长/?=0.2,计算 .y ⑴=1 y(0.2),y(0.4)的近似值，计算过程保留5位小数。 解步长力=0.2,此时/(x,y)=—y—fsiiu 欧拉预报一校正公式为： 预报值兀I = y k + hfg y k) - I J_ 校正值)3=儿+尹(忑，儿)+ fg,儿+1)] 有迭代格式] 预报值儿+] = y k 4-h(-y k -y； sin x k) =y k (0?8-0?2儿sin x k) < h 、—— 2 校止值y如]=儿 +尸[(一片一力sinxJ + LN+i-yl sin.v I+1)] ——?> =儿(°?9一0?1儿sin心)一0?1(儿+| +y；j sin心利) 、"M=0.別=1」)=1 时，Xj=1.2> 有 儿=yo(°?8-O?2yo sinx0) = 1 x (0.8-02x lsin 1) = 0.63171 y(1.2) ?= lx(0.9-0.1xlxsinl)-0.1(0.63171+0.631712sinl.2) = 0.71549 当 T xi=1.2, yi=0.71549 时,x2=1.4,有 y2 =儿(0.8-0?2儿sinXj) = 0.71549x(0.8-02x0.71549sinl.2) =0.47697 y(14) z y2 = 0.71549x(0.9-0.1x0.71549xsin 1.2)-0.1(0.47697+ 0.476972 sin 1.4) =0.52608 V = 8 — 3y 例3写出用四阶龙格一库塔法求解初值问题^ ‘的计算公式，取步长/匸0.2计 b(0) = 2 算y(0.4)的近似值。讣算过程保留4位小数。 解此处.心,刃=8 —3”四阶龙格一库塔法公式为 艰=儿 + % + 2? + 2勺 + ?) 1 h, y n+ y/?A3)： 本例计算公式为: 0 2 呱严儿+三(32?+2?+心 第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 011431313A h A h A h -?=?? ? =?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 令4()f x x =，则 455 1012()5 2 ()(0)()3 h h h h f x dx x dx h A f h A f A f h h ---== -++=? ? 故此时， 101()()(0)()h h f x dx A f h A f A f h --≠-++? 故 101()()(0)()h h f x dx A f h A f A f h --≈-++? 具有3次代数精度。 （2）若 21012()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1014h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 2211163 h h A h A -=+ 从而解得 1143 8383A h A h A h -?=-?? ? =?? ?=?? 令3 ()f x x =，则 22322()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 第一章 绪 论 1．设0x >,x 的相对误差为δ，求ln x 的误差。 解：近似值*x 的相对误差为* **** r e x x e x x δ-=== 而ln x 的误差为()1ln *ln *ln **e x x x e x =-≈ 进而有(ln *)x εδ≈ 2．设x 的相对误差为2%，求n x 的相对误差。 解：设()n f x x =，则函数的条件数为'()||() p xf x C f x = 又1'()n f x nx -= , 1 ||n p x nx C n n -?∴== 又((*))(*)r p r x n C x εε≈? 且(*)r e x 为2 ((*))0.02n r x n ε∴≈ 3．下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位，试指 出它们是几位有效数字：*1 1.1021x =,*20.031x =, *3385.6x =, *456.430x =,*5 7 1.0.x =? 解：*1 1.1021x =是五位有效数字； *20.031x =是二位有效数字； *3385.6x =是四位有效数字； *456.430x =是五位有效数字； *57 1.0.x =?是二位有效数字。 4．利用公式(2.3)求下列各近似值的误差限：(1) *** 124x x x ++,(2) ***123x x x ,(3) **24 /x x . 其中****1234,,,x x x x 均为第3题所给的数。 解： *4 1*3 2*13*3 4*1 51 ()102 1()102 1()102 1()102 1()102x x x x x εεεεε-----=?=?=?=?=? ***124***1244333 (1)() ()()() 111101010222 1.0510x x x x x x εεεε----++=++=?+?+?=? ***123*********123231132143 (2)() ()()() 1111.10210.031100.031385.610 1.1021385.610222 0.215 x x x x x x x x x x x x εεεε---=++=???+???+???≈ **24****24422 *4 33 5 (3)(/)()() 110.0311056.430102256.43056.430 10x x x x x x x εεε---+≈??+??=?= 5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为343V R π= 则何种函数的条件数为 2 3 '4343 p R V R R C V R ππ=== (*)(*)3(*)r p r r V C R R εεε∴≈= 又(*)1r V ε= 数值分析讲义 第三章线性方程组的解法 §3.0 引言 §3.1 雅可比(Jacobi)迭代法 §3.2 高斯-塞德尔(Gauss-Seidel)迭代法 §3.3 超松驰迭代法§3.7 三角分解法 §3.4 迭代法的收敛性§3.8 追赶法 §3.5 高斯消去法§3.9 其它应用 §3.6 高斯主元素消去法§3.10 误差分析 §3 作业讲评3 §3.11 总结 §3.0 引言 重要性：解线性代数方程组的有效方法在计算数学和科学计算中具有特殊的地位和作用.如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题. 分类：线性方程组的解法可分为直接法和迭代法两种方法. (a) 直接法：对于给定的方程组，在没有舍入误差的假设下，能在预定的运算次数内求得精确解.最基本的直接法是Gauss消去法，重要的直接法全都受到Gauss消去法的启发.计算代价高. (b) 迭代法：基于一定的递推格式，产生逼近方程组精确解的近似序列.收敛性是其为迭代法的前提，此外，存在收敛速度与误差估计问题.简单实用,诱人. §3.1 雅可比Jacobi 迭代法 (AX =b ) 1 基本思想: 与解f (x )=0 的不动点迭代相类似,将AX =b 改写为X =BX +f 的形式,建立雅可比方法的迭代格式：X k +1=BX (k )+f ,其中，B 称为迭代矩阵.其计算精度可控,特别适用于求解系数为大型稀疏矩阵(sparse matrices)的方程组. 2 问题: (a) 如何建立迭代格式？ (b) 向量序列{X k }是否收敛以及收敛条件? 3 例题分析: 考虑解方程组??? ??=+--=-+-=--2.453.82102 .72103 21321321x x x x x x x x x (1) 其准确解为X *={1, 1.2, 1.3}. 建立与式(1)相等价的形式： ??? ??++=++=++=84.02.01.083.02.01.072 .02.01.02 13312321x x x x x x x x x (2) 据此建立迭代公式： ?????++=++=++=+++84 .02.01.083.02.01.072.02.01.0)(2)(1)1(3 )(3 )(1)1(23)(2)1(1k k k k k k k k k x x x x x x x x x (3) 取迭代初值0) 0(3 )0(2)0(1===x x x ,迭代结果如下表. JocabiMethodP31.cpp 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27 第一章 绪论（12） 1、设0>x ，x 的相对误差为δ，求x ln 的误差。 [解]设0*>x 为x 的近似值，则有相对误差为δε=)(*x r ，绝对误差为**)(x x δε=，从而x ln 的误差为δδεε=='=* ****1)()(ln )(ln x x x x x ， 相对误差为* * ** ln ln ) (ln )(ln x x x x r δ εε= = 。 2、设x 的相对误差为2%，求n x 的相对误差。 [解]设*x 为x 的近似值，则有相对误差为%2)(*=x r ε，绝对误差为**%2)(x x =ε，从而n x 的误差为n n x x n x n x x n x x x ** 1 *** %2%2) ()()()(ln * ?=='=-=εε， 相对误差为%2) () (ln )(ln *** n x x x n r == εε。 3、下列各数都是经过四舍五入得到的近似数，即误差不超过最后一位的半个单位，试指出它们是几位有效数字： 1021.1*1=x ，031.0*2=x ，6.385*3=x ，430.56*4=x ，0.17*5 ?=x 。 [解]1021.1*1 =x 有5位有效数字；0031.0* 2=x 有2位有效数字；6.385*3=x 有4位有效数字；430.56* 4 =x 有5位有效数字；0.17*5?=x 有2位有效数字。 4、利用公式（3.3）求下列各近似值的误差限，其中* 4*3*2*1,,,x x x x 均为第3题所给 的数。 （1）* 4*2*1x x x ++； [解]3 334* 4*2*11** *4*2*1*1005.1102 1 10211021)()()()()(----=?=?+?+?=++=? ??? ????=++∑x x x x x f x x x e n k k k εεεε； （2）* 3*2 *1x x x ； 一 实验目的 1．掌握求解线性方程组的高斯消元法及列主元素法； 2. 掌握求解线性方程组的克劳特法； 3. 掌握求解线性方程组的平方根法。 二 实验内容 1．用高斯消元法求解方程组（精度要求为610-=ε）： 1231231 233272212240x x x x x x x x x -+=??-+-=-??-+=? 2．用克劳特法求解上述方程组（精度要求为610-=ε）。 3. 用平方根法求解上述方程组（精度要求为610-=ε）。 4. 用列主元素法求解方程组（精度要求为610-=ε）： 1231231 233432222325x x x x x x x x x -+=??-+-=??--=-? 三 实验步骤（算法）与结果 1. 程序代码（Python3.6）： import numpy as np def Gauss(A,b): n=len(b) for i in range(n-1): if A[i,i]!=0: for j in range(i+1,n): m=-A[j,i]/A[i,i] A[j,i:n]=A[j,i:n]+m*A[i,i:n] b[j]=b[j]+m*b[i] for k in range(n-1,-1,-1): b[k]=(b[k]-sum(A[k,(k+1):n]*b[(k+1):n]))/A[k,k] print(b) 运行函数： >>> A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=np.float) >>> b=np.array([7,-1,0],dtype=np.float) >>> x=Gauss(A,b) 输出： 结果：解得原方程的解为x1=3.5,x2=-1,x3=-2.25 2 程序代码（Python3.6）： import numpy as np A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=float) L=np.array([[1,0,0],[0,1,0],[0,0,1]],dtype=float) U=np.array([[0,0,0],[0,0,0],[0,0,0]],dtype=float) b=np.array([7,-1,0],dtype=float) y=np.array([0,0,0],dtype=float) x=np.array([0,0,0],dtype=float) def LU(A): n=len(A[0]) i=0 while i 9．正方形的边长大约为了100cm ，应怎样测量才能使其面积误差不超过2 1cm ？ 解：正方形的面积函数为2 ()A x x = p7 当*100x =时，若(*)1A ε≤, 则21 (*)102 x ε-≤ ? 故测量中边长误差限不超过0.005cm 时，才能使其面积误差不超过2 1cm 第二章 插值法p48 1．当1,1,2 x =-时，()0,3,4f x =-, 分别用单项式基底、拉格朗日基底、牛顿基底求() f x 的二次插值多项式。 解： 0120121200102021101201220211,1,2, ()0,()3,()4;()()1 ()(1)(2)()()2()()1 ()(1)(2) ()()6 ()()1 ()(1)(1) ()()3 x x x f x f x f x x x x x l x x x x x x x x x x x l x x x x x x x x x x x l x x x x x x x ==-===-=--==-+-----==------= =-+-- 则二次拉格朗日插值多项式为 2 20 ()()k k k L x y l x ==∑ 0223()4() 14 (1)(2)(1)(1)23537623 l x l x x x x x x x =-+=---+-+=+- 2．给出()ln f x x =的数值表 用线性插值及二次插值计算的近似值。 解：由表格知， 本科生课程设计报告实习课程数值分析 学院名称管理科学学院 专业名称信息与计算科学 学生姓名 学生学号 指导教师 实验地点 实验成绩 二〇一六年六月二〇一六年六 摘要 常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展. 关键词：数值解法；常微分 1. 实验内容和要求 常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。 要求：分别取步长h = ，，，采用改进的Euler 方法、4阶经典龙格-库塔R －K 方法和4阶Adams 预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。其中多步法需要的初值由经典R-K 法提供。运算时，取足以表示计算精度的有效位。 2. 算法说明 函数及变量说明 表1 函数及变量定义 1、欧拉方法： 1()()(,())i i i i y x y x hf x y x +=+ （i=0，1，2，3，......n-1） (0)y a = (其中a 为初值) 2、改进欧拉方法： ~ 1~111()()(,()) ()()[(,())(,())] 2 (0)i i i i i i i i i i y x y x hf x y x h y x y x f x y x f x y x y a ++++=+=++=（i=0，1，2......n-1） (其中a 为初值) 3、经典K-R 方法： 1 1213 243()6 (,)(,)22(,)22 (,) i i i i i i i i i i h y y K f x y h h K f x y K h h K f x y K K f x h y hK +? =+?? =???=++?? ? =++??=++?? 4、4阶adams 预测-校正方法 预测： 校正： Adsms 内插外插公式联合使用称为Adams 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。计算时需要注意以下两点： 1、外插公式为显式，内插公式为隐式。故用内插外插公式时需要进行迭代。 2、这种预测-校正法是四步法，计算Yn+1时，不但用到前一步信息，而且要用到更前三步信息1-n f ，2-n f 3-n f ，因此它不是自动开始的，实际计算时必须借助某种单步法，用Runge-Kutta 格式为它提供初始值 3 21,,y y y ，依据局部截断误差公式得： 第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h A h -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 01 1431313A h A h A h -?=?? ?=?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++= 故 101()()(0)()h h f x dx A f h A f A f h --=-++? 成立。 令4 ()f x x =，则 第一章 绪论 3．下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位，试指 出它们是几位有效数字：*1 1.1021x =,*20.031x =, *3385.6x =, * 456.430x =,*57 1.0.x =? 解：*1 1.1021x =是五位有效数字； *20.031x =是二位有效数字； *3385.6x =是四位有效数字； *456.430x =是五位有效数字； *57 1.0.x =?是二位有效数字。 4．利用公式(2.3)求下列各近似值的误差限：(1) ***124x x x ++,(2) ***123x x x ,(3) **24/x x . 其中****1234 ,,,x x x x 均为第3题所给的数。 解： *4 1* 3 2* 13* 3 4* 1 51 ()1021()1021()1021()1021()102 x x x x x εεεεε-----=?=?=?=?=? *** 124***1244333 (1)()()()() 1111010102221.0510x x x x x x εεεε----++=++=?+?+?=? *** 123*********123231132143 (2)() ()()() 111 1.10210.031100.031385.610 1.1021385.610222 0.215 x x x x x x x x x x x x εεεε---=++=???+???+???≈ ** 24**** 24422 * 4 33 5 (3)(/) ()() 11 0.0311056.430102256.43056.430 10x x x x x x x εεε---+≈ ??+??= ?= 5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为34 3 V R π= 则何种函数的条件数为 2 3'4343 p R V R R C V R ππ=== (*)(*)3(*)r p r r V C R R εεε∴≈= 又(*)1r V ε= 故度量半径R 时允许的相对误差限为1 (*)10.333 r R ε=?≈ 6．设028Y = ，按递推公式1n n Y Y -= （n=1,2,…） 计算到100Y 27.982≈（5位有效数字），试问计算100Y 将有多大误差？ 解：1n n Y Y -=- 10099Y Y ∴=- 9998Y Y = 9897Y Y =-…… 10Y Y =- 依次代入后，有1000100Y Y =- 即1000Y Y = 27.982, 100027.982Y Y ∴=- ┊┊┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊ 常微分方程的几种数值解法 数学与应用数学肖振华指导教师张秀艳 【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】常微分方程数值解法MA TLAB 误差分析 【Abstract】 Many phenomena in nature and engineering can be attributed to the definite solution of the problem for differential equations. Among them, the ordinary differential equation solving is an important foundation for the content of the differential equations. However, many of the differential equations are often difficult to obtain accurate analytical expression .At this time, the numerical solution provides a good idea. For the Numerical Solution of Ordinary Differential Equations in this article, we focuses on some commonly used numerical solution, such as the Euler method, improved Euler method, Runge-Kutta method, Adams predictor corrector method as well as newer spectral methods. Through specific examples, combined with MATLAB solving and drawing, we initially know the solution process of general numerical solution of ordinary differential equations . At the same time, according to the error analysis of various methods , everyone has an intuitive feel of the characteristics and scope of the various methods. 【Keywords】Ordinary Differential Equations Numerical Solution MATLAB error analysis 《计算机数学基础》数值部分第五单元辅导 14 常微分方程的数值解法 一、重点内容 1. 欧拉公式： ),...,,,() ,()(1-210=???+=+=≈01+1+n k kh x x y x hf y y x y k k k k k k 局部截断误差是O (h 2)。 2. 改进欧拉公式： 预报－校正公式： ?? ???++=+=++++)],(),([2),(1111k k k k k k k k k k y x f y x f h y y y x hf y y 校正值 预报值 即 ))],(,(),([2 11k k k k k k k k y x hf y x f y x f h y y +++ =++ 或表成平均的形式： ??? ? ? ???? +21=+=+=1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y 改进欧拉法的局部截断误差是O (h 3) 3. 龙格－库塔法 二阶龙格－库塔法的局部截断误差是O (h 3) 三阶龙格－库塔法的局部截断误差是O (h 4) 四阶龙格?库塔法公式： )22(6 43211κκκκ++++ =+h y y k k 其中 κ1=f (x k ,y k )；κ2=f (x n +12h ，y k +21h κ1)；κ3=f (x k +12h ，y n +2 1 h κ2)；κ4=f (x k +h ，y k +h κ3) 四阶龙格－库塔法的局部截断误差是O (h 5)。 二、实例 例1 用欧拉法解初值问题???1 =060≤≤0--='2)() .(y x xy y y ，取步长h =0.2。计算过程保留 4位小数。 解h =0.2, f (x )=－y －xy 2。首先建立欧拉迭代格式 ) ,,)((.),(210=-420=--=+=2 1+k y x y y hx hy y y x hf y y k k k k k k k k k k k 当k =0，x 1=0.2时，已知x 0=0,y 0=1，有 实验名称： 实验五 线性方程组的数值解法 指导教师： 数值分析实验组 实验时数： 2 实验设备：安装了Matlab 、C ++、VF 软件的计算机 实验日期：2014年 月 日 实验地点： 第五教学楼北802或902 实验目的： 1. 掌握线性方程组的迭代法和直接法的基本思想和基本步骤。 2. 理解各类数值解法的优缺点，并能自行编程求解。 3. 认识迭代法收敛的含义以及迭代法初值和方程组系数矩阵性质对收敛速度的影响，了解求解病态线性方程组的方法。 实验准备： 1. 在开始本实验之前，请回顾教科书的相关内容； 2. 需要一台准备安装Windows XP Professional 操作系统和装有数学软件的计算机。 实验内容及要求 A 题 考虑方程组b Hx =的求解，其中系数矩阵H 为Hilbert 阵， n j i j i h h H j i n n j i ,,2,1,,1 1 ,)(,, =-+= =? 这是一个著名的病态问题。通过首先给定解（例如取各个分量均为1）再计算出右端的办法给出确定的问题。 （1） 选择问题的维数为6，分别用列主元Gauss 消去法、Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？ （2） 逐步增大问题的维数，仍然用上述方法来解它们，计算的结果如何？计算的结果说明了什么？分析产生结果的原因。 说明：实验过程应包括对问题的简要分析、求解方法、求解步骤、程序及其必要的图表等内容。 实验过程： 本实验所选题为A 题 实验分析：b Hx =，H 矩阵可由Matlab 直接给出，为了设定参考解，首先设x 为分量全 为1的向量，求出b ，然后将H 和b 作为已知量，求x ，与设定的参考解对比。 对于列主元Gauss 消去法，Jacobi 迭代，Gauss-Seidel 迭代，SOR 迭代法，去迭取迭代初值 第一章 绪论 1．设0x >,x 的相对误差为δ，求ln x 的误差。 解：近似值*x 的相对误差为* **** r e x x e x x δ-= = = 而ln x 的误差为()1 ln *ln *ln ** e x x x e x =-≈ 进而有(ln *)x εδ≈ 2．设x 的相对误差为2%，求n x 的相对误差。 解：设()n f x x =，则函数的条件数为'() | |() p xf x C f x = 又1 '()n f x nx -=, 1 ||n p x nx C n n -?∴== 又 ((*))(*)r p r x n C x εε≈? 且(*)r e x 为2 ((*))0.02n r x n ε∴≈ 3．下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位，试指 出它们是几位有效数字：*1 1.1021x =,*20.031x =, *3385.6x =, * 456.430x =,*57 1.0.x =? 解：* 1 1.1021x =是五位有效数字； *20.031x =是二位有效数字； *3385.6x =是四位有效数字； *456.430x =是五位有效数字； *57 1.0.x =?是二位有效数字。 4．利用公式(2.3)求下列各近似值的误差限：(1) ***124x x x ++,(2) ***123x x x ,(3) **24/x x . 其中**** 1234,,,x x x x 均为第3题所给的数。 解： *4 1* 3 2* 13* 3 4* 1 51()1021()1021()1021()1021()102 x x x x x εεεεε-----=?=?=?=?=? *** 124***1244333 (1)()()()() 1111010102221.0510x x x x x x εεεε----++=++=?+?+?=? *** 123*********123231132143 (2)() ()()() 111 1.10210.031100.031385.610 1.1021385.610222 0.215 x x x x x x x x x x x x εεεε---=++=???+???+???≈ ** 24**** 24422 *4 33 5 (3)(/) ()() 11 0.0311056.430102256.43056.430 10x x x x x x x εεε---+≈ ??+??= ?= 5计算球体积要使相对误差限为1，问度量半径R 时允许的相对误差限是多少？ 解：球体体积为343 V R π= 则何种函数的条件数为 2 3'4343 p R V R R C V R ππ=== (*)(*)3(*)r p r r V C R R εεε∴≈= 又 (*)1r V ε=%1数值分析常微分方程的数值解法
数值分析第四章数值积分与数值微分习题答案
常微分方程数值解法
数值分析第五版全答案chap1
数值分析讲义线性方程组的解法
常微分方程的数值解
最新数值分析课程第五版课后习题答案(李庆扬等)1
数值分析实验2_求解线性方程组直接法
数值分析常微分数值解的求法C语言
数值分析第五版全答案
数值分析第五版答案
常微分方程的几种数值解法
数值分析常微分方程的数值解法
数值分析— 线性方程组的数值解法
数值分析第五版答案(全)
相关主题
文本预览