(完整版)数值计算方法上机实习题答案
计算方法上机实习题
数值计算方法上机实习题1. 设⎰+=105dx xx I nn , (1) 由递推公式nI I n n 151+-=-,从0=0.1822I , 0=0.1823I 出发,计算20I ; (2) 20=0I ,20=10000I , 用nI I n n 51511+-=-,计算0I ;(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
解:(1)程序如下: clear all clc I=0.1822; %题中的已知数据 for n=1:20; I=(-5)*I+1/n; %由递推公式所得 end fprintf('I20=%f\n',I) M=0.1823; %与I 的计算结果形成对比 for i=1:20; M=(-5)*M+1/i; %由递推公式所得 end fprintf('M20=%f\n',M) 输出结果为: I20=-11592559237.912731 M20=-2055816073.851284 (2)程序如下: clear all clc I=0; %赋予I20的初始值 for n=0:19; I=(-1/5)*I+1/(5*(20-n)); %有递推公式得 end fprintf('I0=%f\n',I) M=10000; for i=0:19; M=(-1/5)*M+1/(5*(20-i));%有递推公式得 end fprintf('M0=%f\n',M) 输出结果为: I0=0.182322 M0=0.182322(3)由输出结果可看出第一种算法为不稳定算法,第二中算法为稳定算法。
由于误差*000***21111120115(5)5()555nn n n n n n n n n e I I e I I I I I I e e e n n------=-=-=-+--+=-===第一种算法为正向迭代算法,每计算一步误差增长5倍,虽然所给的初始值很接近,随着n 的增大,误差也越来越大。
(完整版)数值计算方法上机实习题答案
(完整版)数值计算⽅法上机实习题答案1.设?+=105dx xx I nn ,(1)由递推公式nI I n n 151+-=-,从0I 的⼏个近似值出发,计算20I ;解:易得:0I =ln6-ln5=0.1823, 程序为:I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为:20I = -3.0666e+010 (2)粗糙估计20I ,⽤nI I n n 515111+-=--,计算0I ;因为 0095.056 0079.01020201020≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为:I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083(3)分析结果的可靠性及产⽣此现象的原因(重点分析原因)。
⾸先分析两种递推式的误差;设第⼀递推式中开始时的误差为000I I E '-=,递推过程的舍⼊误差不计。
并记nn n I I E '-=,则有01)5(5E E E n n n -==-=-Λ。
因为=20E 20020)5(I E >>-,所此递推式不可靠。
⽽在第⼆种递推式中n n E E E )51(5110-==-=Λ,误差在缩⼩,所以此递推式是可靠的。
出现以上运⾏结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。
2.求⽅程0210=-+x e x的近似根,要求41105-+?<-k k x x ,并⽐较计算量。
(1)在[0,1]上⽤⼆分法;程序:a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果:c =0.0903(2)取初值00=x ,并⽤迭代1021x k e x -=+;程序:x=0; a=1;while abs(x-a)>5*1e-4 a=x;x=(2-exp(x))/10; end x结果:x =0.0905(3)加速迭代的结果;程序:x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果:x =0.0995(4)取初值00=x ,并⽤⽜顿迭代法;程序:x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x结果: x =0.0905(5)分析绝对误差。
计算方法与实习答案1-2
绪论
习题1——10:设 f ( x) = 8 x 5 − 0.4 x 4 + 4 x 3 − 9 x + 1 用秦九韶法求f(3)。 解:
8 − 0.4
24 8 23.6
0
−9
1
x=3
70.8 74.8
224.4 224.4
673.2 664.2
1992.6 1993.6
∴ f(3)=1993.6
第一章 绪论 练习
1.《计算方法》课程主要研究以计算 机为工具的 数值 分析方法 ,并评价 该算法的计算误差。 2.近似值作四则运算后的绝对误差限 公式为 ε ( x1 − x2 ) ≤ ε ( x1 ) + ε ( x2 ) ,近似值 1.0341的相对误差限不大于 1 ×10−2 , 则它至少有三位有效数字。 4
ln(103 ) ∴k ≥ ln(2) ≥ 9.965
2 2 2
∴需二分10次 需二分 次
方程求根——二分法
习题2——2:用二分法求方程2e-x-sinx=0在区 间[0,1]内的1个实根,要求3位有效数字。
解:1)判断是否在该区间有且仅有一个根 f(0)=2>0,f(1)=2/e-sin1≈-0.1<0, f’(x)=-2e-x-cosx,f’=-3,-2/e-cos1<0 2)判断二分次数 由(b-a)/2k+1=1/2k+1≤1/2*10-3,解得k≥3ln10/ln2≥9.965, 所以需要二分10次,才能满足精度要求。
∴ x≈2.981
方程求根
f (xk )(xk − xk −1) xk +1 = xk − f (xk ) − f (xk −1)
习题2——11:用割线法求方程x3-2x-5=0的根,要 求精确到4位有效数字,取x0=2, x1=2.2。
(完整版)数值计算方法试题及答案
数值计算方法试题一一、 填空题(每空1分,共17分)1、如果用二分法求方程043=-+x x 在区间]2,1[内的根精确到三位小数,需对分( )次。
2、迭代格式)2(21-+=+k k k x x x α局部收敛的充分条件是α取值在( )。
3、已知⎪⎩⎪⎨⎧≤≤+-+-+-≤≤=31)1()1()1(2110)(233x c x b x a x x x x S 是三次样条函数,则a =( ),b =( ),c =( )。
4、)(,),(),(10x l x l x l n 是以整数点n x x x ,,,10 为节点的Lagrange 插值基函数,则∑==nk kx l0)(( ),∑==nk k jk x lx 0)(( ),当2≥n 时=++∑=)()3(204x l x xk k n k k( )。
5、设1326)(247+++=x x x x f 和节点,,2,1,0,2/ ==k k x k 则=],,,[10n x x x f 和=∆07f。
6、5个节点的牛顿-柯特斯求积公式的代数精度为 ,5个节点的求积公式最高代数精度为 。
7、{}∞=0)(k kx ϕ是区间]1,0[上权函数x x =)(ρ的最高项系数为1的正交多项式族,其中1)(0=x ϕ,则⎰=14)(dx x x ϕ 。
8、给定方程组⎩⎨⎧=+-=-221121b x ax b ax x ,a 为实数,当a 满足 ,且20<<ω时,SOR 迭代法收敛。
9、解初值问题00(,)()y f x y y x y '=⎧⎨=⎩的改进欧拉法⎪⎩⎪⎨⎧++=+=++++)],(),([2),(]0[111]0[1n n n n n n n n n n y x f y x f h y y y x hf y y 是阶方法。
10、设⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=11001a a a a A ,当∈a ( )时,必有分解式T LL A =,其中L 为下三角阵,当其对角线元素)3,2,1(=i l ii 满足( )条件时,这种分解是唯一的。
上机数值计算练习题及答案.docx
习题31、在MATLAB 中,先运行指令A=magic(3), B=[l,2,l;3,4,3;5,6,7], C=reshape(l:6,3,2)生成阵列A 珂,B3X2,C3X2 ,然后根据运行结果回答以下问题:(1)计算A*B,B*A,这两个乘积相同吗?(2)计算A\B, B/A,左除、右除结果相同吗?(3)计算B(:,[1,2]).*C和C.*B(:, [1,2]),这两个乘积相同吗?(4)计算A\A和AAA,这两个计算结果相同吗?(5)计算A\eye(3)和inv(A),这两个计算结果相同吗?(提示:根据对计算结果的目测回答问题)解答如下:运行指令:A=magic(3), B=[l,2,l;3,4,3;5,6,7], C=reshape(l:6,3,2)得到结果:8 1 63 5 74 9 2B =1 2 13 4 35 6 7C =1 42 53 6(1)计算A*B,并得到结果如下:» A*Bans =41 56 5353 68 6741 56 45计算B*A, 并得到结果如下:»B*Aans =18 20 2248 50 5286 98 86从以上计算结果可以得出结论:A*BJ (2)计算A\B ,并得到结果如下:» A\Bans =0.0333 0.1000 0.16110.5333 0.6000 0.74440.0333 0.1000 -0.1722计算B/A, 并得到结果如下:B/Aans =0.0056 0.0889 0.17220.1389 0.2222 0.30560.2333 0.7333 0.2333 与B*A这两个乘积不同。
从以上计算结杲可以得出结论:左除、右除结杲不同。
(3)计算B(:,[1,2]).弋,并得到结果如下:A =» B(:,[1,2]).*C ans =1 8 6 20 15 36计算C.*B(:, [1,2]),并得到结果如下: » CFB(:, [1,2]) ans =1 6 20 15 36从以上计算结果可以得出结论:B(: J1,2]).*C 和C ・*B(:, [1,2])的两个乘积相同。
数值分析上机实习题及答案.docx
方詡文金兴:爭[数值分析]2017-2018第二学期上机实习题1:编程计算亍丄,其中C= 4. 4942x10307,给出并观察计算结心C"果,若有问题,分析之。
解:mat lab 编程如下:E) funct ion diy i ti formatlong g;n 二input C 输入ii 值= c= 4.4942E307; sum 二0; s 二 0;E3 for i = l:n s = l/ (c*i);>> diyiti 输入n 值:10 104.6356e-308 >> diyiti输入ri 值:1001004.6356e-308 >> diyiti 输入n 值:1000 10004.6356e-308 >> diyiti揄入n 值* 1000001000004・ 6356e-308 >> diyiti输入n 值;1000000001000000004.6356e-308图二:运行结果Mat lab 中,forma t long g 对双精度,显示15位定点或浮点格式,由上图 可知,当输入较小的n 值5分别取10, 100, 1000, 100000, 100000000)的时候, 结果后面的指数中总是含有e-308,这和题目中的C 值很相似,我认为是由于分 母中的C 值相对于n 值过大,出现了 “大数吃小数”的现彖,这是不符合算法原 则的。
2:利用牛顿法求方程-1^ = 2于区间241的根,考虑不同初值下牛顿法的收敛情况。
解:牛顿法公式为:利用mat lab 编程function di2ti21 3i=l ;2 2.85208156699784 xO 二input ('输入初值x0:‘ );A 二[i x0];3 2.55030468822809 t=x0+ (x0-log (xO) -2) /(1-1/xO) ; %迭代函数4 1.91547247100476 三 while (abs (t _x0)>0.01)i=i+l; 5 0.37867158538991 xO 二 t; 6 0.774964959780275 A = [A;i xO];t =x0+(x0-log(xO)-2)/(1-1/xO): 7 4.11574081641933 cnd| 8 5.04162436446126 disp (A);96.81782826645596当输入初值二3的时候并不能收敛。
数值计算方法试题和答案解析
数值计算方法试题一一、 填空题(每空1分,共17分) 1、如果用二分法求方程043=-+x x 在区间]2,1[内的根精确到三位小数,需对分( )次。
2、迭代格式)2(21-+=+k k k x x x α局部收敛的充分条件就是α取值在( )。
3、已知⎪⎩⎪⎨⎧≤≤+-+-+-≤≤=31)1()1()1(2110)(233x c x b x a x x x x S 就是三次样条函数,则a =( ),b =( ),c =( )。
4、)(,),(),(10x l x l x l n Λ就是以整数点n x x x ,,,10Λ为节点的Lagrange 插值基函数,则∑==nk kx l0)(( ),∑==nk k jk x lx 0)(( ),当2≥n 时=++∑=)()3(204x l x xk k n k k( )。
5、设1326)(247+++=x x x x f 与节点,,2,1,0,2/Λ==k k x k 则=],,,[10n x x x f Λ 与=∆07f。
6、5个节点的牛顿-柯特斯求积公式的代数精度为 ,5个节点的求积公式最高代数精度为 。
7、{}∞=0)(k kx ϕ就是区间]1,0[上权函数x x =)(ρ的最高项系数为1的正交多项式族,其中1)(0=x ϕ,则⎰=14)(dx x x ϕ 。
8、给定方程组⎩⎨⎧=+-=-221121b x ax b ax x ,a 为实数,当a 满足 ,且20<<ω时,SOR 迭代法收敛。
9、解初值问题00(,)()y f x y y x y '=⎧⎨=⎩的改进欧拉法⎪⎩⎪⎨⎧++=+=++++)],(),([2),(]0[111]0[1n n n n n n n n n n y x f y x f h y y y x hf y y 就是阶方法。
10、设⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=11001a a a a A ,当∈a ( )时,必有分解式T LL A =,其中L为下三角阵,当其对角线元素)3,2,1(=i l ii 满足( )条件时,这种分解就是唯一的。
数值计算方法答案
数值计算方法答案(总71页)--本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--数值计算方法习题一(2)习题二(6)习题三(15)习题四(29)习题五(37)习题六(62)习题七(70)2009.9,9习题一1.设x >0相对误差为2%4x 的相对误差。
解:由自变量的误差对函数值引起误差的公式:(())(())'()()()()f x xf x f x x f x f x δδ∆=≈得(1)()f x =11()()*2%1%22x x δδδ≈===;(2)4()f x x =时444()()'()4()4*2%8%x x x x x xδδδ≈===2.设下面各数都是经过四舍五入得到的近似数,即误差不超过最后一位的半个单位,试指出他们各有几位有效数字。
(1)12.1x =;(2)12.10x =;(3)12.100x =。
解:由教材9P 关于1212.m nx a a a bb b =±型数的有效数字的结论,易得上面三个数的有效数字位数分别为:3,4,53.用十进制四位浮点数计算 (1)++; (2)+(+)哪个较精确 解:(1)++ ≈21((0.3197100.245610)0.1352)fl fl ⨯+⨯+=2(0.3443100.1352)fl ⨯+ =210⨯ (2)+(+) 21(0.319710(0.245610))fl fl ≈⨯+⨯= 21(0.3197100.259110)fl ⨯+⨯=210⨯易见++=210⨯,故(2)的计算结果较精确。
4.计算正方形面积时,若要求面积的允许相对误差为1%,测量边长所允许的相对误差限为多少解:设该正方形的边长为x ,面积为2()f x x =,由(())(())'()()()()f x xf x f x x f x f x δδ∆=≈解得(())()()'()f x f x x xf x δδ≈=2(())(())22f x x f x x xδδ==%5.下面计算y 的公式哪个算得准确些为什么(1)已知1x <<,(A )11121xy x x -=-++,(B )22(12)(1)x y x x =++; (2)已知1x >>,(A )y=,(B )y = (3)已知1x <<,(A )22sin x y x=,(B )1cos 2xy x -=;(4)(A)9y =B )y =解:当两个同(异)号相近数相减(加)时,相对误差可能很大,会严重丧失有效数字;当两个数相乘(除)时,大因子(小除数)可能使积(商)的绝对值误差增大许多。
数值计算课后习题答案(全)
习 题 一 解 答1.取3.14,3.15,227,355113作为π的近似值,求各自的绝对误差,相对误差和有效数字的位数。
分析:求绝对误差的方法是按定义直接计算。
求相对误差的一般方法是先求出绝对误差再按定义式计算。
注意,不应先求相对误差再求绝对误差。
有效数字位数可以根据定义来求,即先由绝对误差确定近似数的绝对误差不超过那一位的半个单位,再确定有效数的末位是哪一位,进一步确定有效数字和有效数位。
有了定理2后,可以根据定理2更规范地解答。
根据定理2,首先要将数值转化为科学记数形式,然后解答。
解:(1)绝对误差:e(x)=π-3.14=3.14159265…-3.14=0.00159…≈0.0016。
相对误差:3()0.0016()0.51103.14r e x e x x -==≈⨯有效数字:因为π=3.14159265…=0.314159265…×10,3.14=0.314×10,m=1。
而π-3.14=3.14159265…-3.14=0.00159…所以│π-3.14│=0.00159…≤0.005=0.5×10-2=21311101022--⨯=⨯所以,3.14作为π的近似值有3个有效数字。
(2)绝对误差:e(x)=π-3.15=3.14159265…-3.14=-0.008407…≈-0.0085。
相对误差:2()0.0085()0.27103.15r e x e x x --==≈-⨯有效数字:因为π=3.14159265…=0.314159265…×10,3.15=0.315×10,m=1。
而π-3.15=3.14159265…-3.15=-0.008407…所以│π-3.15│=0.008407……≤0.05=0.5×10-1=11211101022--⨯=⨯所以,3.15作为π的近似值有2个有效数字。
(3)绝对误差:22() 3.14159265 3.1428571430.0012644930.00137e x π=-=-=-≈-相对误差:3()0.0013()0.4110227r e x e x x--==≈-⨯有效数字:因为π=3.14159265…=0.314159265…×10, 223.1428571430.3142857143107==⨯,m=1。
数值分析计算实习题答案
数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。
在数值分析的学习过程中,实习题是一种重要的学习方式,通过实践来巩固理论知识,并培养解决实际问题的能力。
本文将为大家提供一些数值分析计算实习题的答案,希望能够帮助大家更好地理解和掌握数值分析的相关知识。
一、插值与拟合1. 已知一组数据点,要求通过这些数据点构造一个一次插值多项式,并求出在某一特定点的函数值。
答案:首先,我们可以根据给定的数据点构造一个一次插值多项式。
假设给定的数据点为(x0, y0), (x1, y1),我们可以构造一个一次多项式p(x) = a0 + a1x,其中a0和a1为待定系数。
根据插值条件,我们有p(x0) = y0,p(x1) = y1。
将这两个条件代入多项式中,可以得到一个方程组,通过求解这个方程组,我们就可以确定a0和a1的值。
最后,将求得的多项式代入到某一特定点,就可以得到该点的函数值。
2. 已知一组数据点,要求通过这些数据点进行最小二乘拟合,并求出拟合曲线的表达式。
答案:最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。
假设给定的数据点为(x0, y0), (x1, y1),我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n,其中a0, a1, ..., an为待定系数。
根据最小二乘拟合原理,我们需要最小化误差平方和E = Σ(yi - f(xi))^2,其中yi为实际数据点的y值,f(xi)为拟合曲线在xi处的函数值。
通过求解这个最小化问题,我们就可以确定拟合曲线的表达式。
二、数值积分1. 已知一个函数的表达式,要求通过数值积分的方法计算函数在某一区间上的定积分值。
答案:数值积分是一种通过将定积分转化为数值求和来近似计算的方法。
假设给定的函数表达式为f(x),我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi),其中wi为权重系数,xi为待定节点。
计算方法与实习上机题答案
实习题11用两种不容的顺序计算644834.1100012≈∑=-n n,分析误差的变化(1)顺序计算 源代码:#include<stdio.h> #include<math.h> void main() { double sum=0; int n=1; while(1) { sum=sum+(1/pow(n,2)); if(n%1000==0)printf("sun[%d]=%-30f",n,sum); if(n>=10000)break; n++; } printf("sum[%d]=%f\n",n,sum); }结果:(2)逆序计算 源代码:#include<stdio.h> #include<math.h> void main() { double sum=0; int n=10000; while(1) { sum=sum+(1/pow(n,2));if(n%1000==0) printf("sum[%d]=%-30f",n,sum); if(n<=1)break; n --;} printf("sum[%d]=%f\n",n,sum);}结果:2已知连分数))//(.../(322101n n b a a b a b a b f ++++=利用下面的方法计算f:11)0,...,2,1(,d f n n i d a b d b d i i i i n n =--=+==++写一个程序,读入n,n n b a ,,计算并打印f 源代码:#include<stdio.h> #include<math.h> void main() { int i=0,n; float a[1024],b[1024],d[1024]; printf("please input n,n="); scanf("%d",&n); printf("\nplease input a[1] to a[n]:\n"); for(i=1;i<=n;i++) { printf("a[%d]=",i); scanf("%f",&a[i]);} printf("\nplease input b[0] to b[n]:\n"); for(i=0;i<=n;i++) { printf("b[%d]=",i); scanf("%f",&b[i]); } d[n]=b[n]; for(i=n -1;i>=0;i --) d[i]=b[i]+a[i+1]/d[i+1]; printf("\nf=%f\n",d[0]); } 结果:3给出一个有效的算法和一个无效的算法计算积分⎰=+=10)10,...1,0(14n dx x x y n n源代码:#include<stdio.h> #include<math.h> main() { double y_0=(1/4.0)*log(5),y_1; double y_2=(1.0/55.0+1.0/11.0)/2,y_3; int n=1,m=10; printf("有效算法输出结果:\n"); printf("y[0]=%-20f",y_0);while(1) { y_1=1.0/(4*n)+y_0/(-4.0); printf("y[%d]=%-20f",n,y_1); if(n>=10) break; y_0=y_1; n++; if(n%3==0) printf("\n"); }printf("\n 无效算法的输出结果:\n"); printf("y[10]=%-20f",y_2); while(1) { y_3=1.0/n -4.0*y_2; printf("y[%d]=%-20f",m -1,y_3); if(m<=1) break; y_2=y_3; m --; if(m%2==0) printf("\n"); }}结果:4设∑=-=Nj N j S 2211,已知其精确值为)11123(21+--N N (1)编制按从小到大顺序计算N S 的程序 (2)编制按从小达到的顺序计算N S 的程序(3)按两种顺序分别计算30000100001000,,S S S ,并指出有效位数源代码:#include<stdio.h> main() { int N; double SN[30000]; SN[30000]=(3.0/2.0-1.0/30000.0-1/30001.0)/2.0; for(N=30000;N>=2;N --) SN[N -1]=SN[N]-1.0/(N*N -1); printf("从大到小顺序计算:\nSN[1000]=%f\nSN[10000]=%f\nSN[30000]=%f\n",SN[1000],SN[10000],SN[30000]); SN[2]=(3.0/2-1.0/2.0-1/3.0)/2.0; for(N=3;N<=30000;N++) SN[N]=SN[N -1]+1.0/(N*N -1); printf("从小到大顺序计算:\nSN[1000]=%f\nSN[10000]=%f\nSN[30000]=%f\n",SN[1000],SN[10000],SN[30000]); } 结果:实习题2 1.用牛顿法求下列方程的根2=-x e x01=-x xe 02lg =-+x x源代码:#include <stdio.h> #include <math.h>typedef float (*p)(float ); float ff1(float x) { return x*x -exp(x); }float ff2(float x) { return x*exp(x)-1; }float ff3(float x) { return log(x)+x -2; }float answer(float(*p)(float)) { int k=2; float m=1,n=-1,x2,a,b,c; if (p==ff3)n=2; printf("x[0] = %.4f, x[1] = %.4f, ",m,n); while (1) { if (fabs(m -n)<1e -4) break; a=p(n)*(n -m); b=p(n)-p(m); c=a/b; x2=n -c; m = n; n = x2; printf("x[%d] = %.4f, ",k,x2); k++; if (k%3==0) printf("\n"); }if (k%3!=0) printf("\n");printf("iteration times: %d, roots: %.4f\n ",k-2,n);return 0;}main(){printf("x*x-exp(x),\n");answer(ff1);printf("x*exp(x)-1,\n");answer(ff2);printf("lg(x)+x-2,\n");answer(ff3);return 0;}结果:2.编写一个割线法的程序,求解上述各方程源代码:#include<stdio.h>#include<math.h>float gexian(float,float);float f(float);main(){int i,j;float x1=2.2;float x2=2,x3;scanf("%d",&i);if(i==1)printf("%f",x1); else if(i==2) printf("%f",x2); else { for(j=3;j<=i;j++) { x3=gexian(x1,x2); x1=x2; x2=x3; } printf("%f",gexian(x1,x2)); } }float f(float x) { return (x*x -exp(x)); }float gexian(float x1,float x2) { return (x2-(f(x2)/(f(x2)-f(x1)))*(x2-x1)); }结果:实习题3 1用列主元消去法解下列方程组;⎪⎪⎩⎪⎪⎨⎧=++=-++--=+--=--+43443233312)1(421432143214321x x x x x x x x x x x x x x x ⎪⎪⎩⎪⎪⎨⎧=++--=++-=-+--=-+-4341220332282)2(432132143214321x x x x x x x x x x x x x x x 源程序:#include<stdio.h>#include<math.h>void ColPivot(float*,int,float[]);void ColPivot(float*c,int n,float x[]){int i,j,t,k;float p;for(i=0;i<=n-2;i++){k=i;for(j=i+1;j<=n-1;j++)if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))k=j;if(k!=i)for(j=i;j<=n;j++){p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;}for(j=i+1;j<=n-1;j++){p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));}}for(i=n-1;i>=0;i--){for(j=n-1;j>=i+1;j--)(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));}}void main(){int i;float x[4];float c[4][5]={1,1,0,3,4,2,1,-1,1,1,3,-1,-1,3,-3,-1,2,3,-1,4};ColPivot(c[0],4,x);for(i=0;i<=3;i++)printf("x[%d]=%f\n",i,x[i]);}结果:第(1)题第(2)题2、源代码:#include<stdio.h>void main(){float x[4];int i;float a[4][5]={48,-24,0,-12,4,-24,24,12,12,4,0,6,20,2,-2,-6,6,2,16,-2};void DirectLU(float*,int,float[]);DirectLU(a[0],4,x);for(i=0;i<=3;i++)printf("x[%d]=%f\n",i,x[i]);}void DirectLU(float*u,int n,float x[]){int i,r,k;for(r=0;r<=n-1;r++){for(i=r;r<=n;i++)for(k=0;k<=r-1;k++)*(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i));for(i=r+1;i<=n-1;i++){for(k=0;k<=r-1;k++)*(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+r));*(u+i*(n+1)+r)/=*(u+r*(n+1)+r);}}for(i=n-1;i>=0;i--){for(r=n-1;r>=i+1;r--)*(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*x[r];x[i]=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i));}}实习题41、源代码:#include<stdio.h>float Lagrange(float x[],float y[],float xx,int n) //n为(n+1)次插值;{int i,j;float *a,yy=0;a=new float[n];for(i=0;i<=n-1;i++){a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}delete a;return yy;}void main(){float x[5]={-3.0,-1.0,1.0,2.0,3.0};float y[5]={1.0,1.5,2.0,2.0,1.0};float xx1=-2,xx2=0,xx3=2.75,yy1,yy2,yy3;yy1=Lagrange(x,y,xx1,3);yy2=Lagrange(x,y,xx2,3);yy3=Lagrange(x,y,xx3,3);printf("x1=%-20f,y1=%f\n",xx1,yy1);printf("x2=%-20f,y2=%f\n",xx2,yy2);printf("x3=%-20f,y3=%f\n",xx3,yy3);}结果:2、源代码:#include<stdio.h>float Lagrange(float x[],float y[],float xx,int n) //n为(n+1)次插值;{int i,j;float *a,yy=0;a=new float[n];for(i=0;i<=n-1;i++){a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}delete a;return yy;}void main(){float x[6]={0.30,0.42,0.50,0.58,0.66,0.72};float y[6]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};float xx1=0.46,xx2=0.55,xx3=0.60,yy1,yy2,yy3;yy1=Lagrange(x,y,xx1,6);yy2=Lagrange(x,y,xx2,6);yy3=Lagrange(x,y,xx3,6);printf("x1=%-20f,y1=%f\n",xx1,yy1);printf("x2=%-20f,y2=%f\n",xx2,yy2);printf("x3=%-20f,y3=%f\n",xx3,yy3);}结果:源代码:#include<stdio.h>#define N 3void Difference(float y[],float f[4][4],int n){int k,i;f[0][0]=y[0];f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];for(k=1;k<=n;k++)for(i=0;i<=(N-k);i++)f[i][k]=f[i+1][k-1]-f[i][k-1];return;}void main(){int i,k=1;float a,b=1,m=21.4,t=1.4,f[4][4]={0};float x[5]={20,21,22,23,24};float y[5]={1.30103,1.32222,1.34242,1.36173,1.38021};Difference(y,f,N);a=f[0][0];for(i=1;i<=N;i++){k=k*i;b=b*(t-i+1);a=a+b*f[0][i]/k;}printf("x(k)\n");for (i=0;i<=4;i++)printf( "%-20f",x[i]);printf("\ny(k)\n");for (i=0;i<=4;i++)printf("%-20f",y[i]);for(k=1;k<=3;k++){printf("\nF(%d)\n ",k);for(i=0;i<=(3-k);i++){printf("%-20f",f[i][k]);}}printf ("\n");printf("f(%f)=%-20f",m,a);printf ("\n");结果:实习题52、源代码:#include<stdio.h>#include<math.h>void main(){int i,n;float a[2];float x[15]={1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8},z[15];floaty[15]={33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.50,267.55,280.50,296.65,301.40,310 .40,318.15,325.15};for(n=0;n<=14;n++) //增加了数组z;{z[n]=log(y[n]/x[n]);}void Approx(float[],float[],int,int,float[]);Approx(x,z,15,1,a); //变成一次拟合;//for(i=0;i<=1;i++)//printf("a[%d]=%f\n",i,a[i]);printf("a=exp(a[0])=%f\n",exp(a[0]));printf("b=-a[1]=%f\n",-a[1]); }void Approx(float x[],float y[],int m,int n,float a[]){int i,j,t;float *c=new float[(n+1)*(n+2)];float power(int,float);void ColPivot(float *,int,float[]);for(i=0;i<=n;i++){for(j=0;j<=n;j++){*(c+i*(n+2)+j)=0;for(t=0;t<=m-1;t++)*(c+i*(n+2)+j)+=power(i+j,x[t]);}*(c+i*(n+2)+n+1)=0;for(j=0;j<=m-1;j++)*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);}ColPivot(c,n+1,a);delete c;}void ColPivot(float *c,int n,float x[]){int i,j,t,k;float p;for(i=0;i<=n-2;i++){k=i;for(j=i+1;j<=n-1;j++)if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))k=j;if(k!=i)for(j=i;j<=n;j++){p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;}for(j=i+1;j<=n-1;j++){p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));}}for(i=n-1;i>=0;i--){for(j=n-1;j>=i+1;j--)(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));}float power(int i,float v){float a=1;while(i--)a*=v;return a;}结果:实习题61、源代码:(1)#include<stdio.h>#include<math.h>float Cotes(float(*f)(float),float a,float b,int n){int k;float c,c1=0,c2,c3,c4;float h=(b-a)/n;c2=(*f)(a+h/4);c3=(*f)(a+h/2);c4=(*f)(a+3*h/4);for(k=1;k<=n-1;k++){c1+=(*f)(a+k*h);c2+=(*f)(a+k*h+h/4);c3+=(*f)(a+k*h+h/2);c4+=(*f)(a+k*h+3*h/4);}c=h/90*(7*((*f)(a)+(*f)(b))+14*c1+32*c2+12*c3+32*c4);return c;}float f(float x){return 1/sqrt(1+x*x*x);}void main()int i,n=4;float c;for(i=0;i<=4;i++){c=Cotes(f,0,1,n);printf("C(%d)=%f\n",n,c);n*=2;}}(2)#include<stdio.h>#include<math.h>float Cotes(float(*f)(float),float a,float b,int n){int k;float c,c1=0,c2,c3,c4;float h=(b-a)/n;c2=(*f)(a+h/4);c3=(*f)(a+h/2);c4=(*f)(a+3*h/4);for(k=1;k<=n-1;k++){c1+=(*f)(a+k*h);c2+=(*f)(a+k*h+h/4);c3+=(*f)(a+k*h+h/2);c4+=(*f)(a+k*h+3*h/4);}c=h/90*(7*((*f)(a)+(*f)(b))+14*c1+32*c2+12*c3+32*c4);return c;}float f(float x){ // return 1/sqrt(1+x*x*x);if (x==0)return 1;else return sin(x)/x;}void main(){int i,n=4;float c;for(i=0;i<=4;i++){// c=Cotes(f,0,1,n);c=Cotes(f,0,5,n);printf("C(%d)=%f\n",n,c);n*=2;}}结果:(1)(2)实习题7一、改进欧拉法1、#include<stdio.h>#include<iostream>double Adams (double (*f)(double x, double y),double x0,double y0,double h,int N) {for(int n=0; n<N; ++n) {double x1=x0+h;double yp=y0+h*f(x0,y0);double y1=y0+h*f(x1,yp);y1=(yp+y1)/2.0;printf("ty=%f\n",y1);x0=x1;y0=y1;}}int main(void){double f(double x, double y); double x0=0,y0=0;double x,y,step;long i;step=0.1;Adams(f,x0,y0,0.1,10);}double f(double x, double y) {double r;r=x*x+y*y;return(r);}2、int main(void){double f(double x, double y); double x0=0,y0=0;double x,y,step;long i;step=0.1;Adams(f,x0,y0,0.1,10);}double f(double x, double y) {double r;r=1/(1+y*y);return(r);}3、int main(void){double f(double x, double y); double x0=0,y0=1;double x,y,step;long i;step=0.1;Adams(f,x0,y0,0.1,10);}double f(double x, double y) {double r;r=y-2*x/y;}四阶龙格库塔法1、#include<iostream>#include<stdio.h>using namespace std;//f为函数的入口地址,x0、y0为初值,xn为所求点,step为计算次数double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, long step ){double k1,k2,k3,k4,result;double h=(xn-x0)/step;if(step<=0)return(y0);if(step==1){k1=f(x0,y0);k2=f(x0+h/2, y0+h*k1/2);k3=f(x0+h/2, y0+h*k2/2);k4=f(x0+h, y0+h*k3);result=y0+h*(k1+2*k2+2*k3+k4)/6;}else{double x1,y1;y1=Runge_Kuta(f, x0, y0, xn-h,step-1);k1=f(x1,y1);k2=f(x1+h/2, y1+h*k1/2);k3=f(x1+h/2, y1+h*k2/2);k4=f(x1+h, y1+h*k3);result=y1+h*(k1+2*k2+2*k3+k4)/6;}printf("%lg\n",result);return(result);}int main(void){double f(double x, double y);double x0=0,y0=0,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl; //}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=x*x+y*y;return(r);}2、int main(void){double f(double x, double y);double x0=0,y0=0,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl; //}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=1/(1+y*y);return(r);}3、int main(void){double f(double x, double y);double x0=0,y0=0,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl; //}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=1/(1+y*y);return(r);}二、int main(void){double f(double x, double y);double x0=0,y0=1,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl;//}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=y-2*x/y;return(r);}三、1、void Runge_Kutta(float(*f)(float x,float y),float a,float b,float y0,int N,float yy[]) {float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;for(i=1;i<=3;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/2,y+h*K3);y=y+h*(K1+2*K2+2*K2+2*K3+K4)/6;x=a+i*h;yy[i-1]=y;}}void Adams (float a,float b,int N,float(*f)(float x,float y),float y0){int i;float y1,y2,y,yp,yc,yy[3],h,x;printf("x[0]=%f\ty[0]=%f\n",a,y0);Runge_Kutta(f,a,b,y0,N,yy);y1=yy[0];y2=yy[1];y=yy[2];h=(b-a)/N;for(i=1;i<=3;i++)printf("x[%d]=%f\ty[%d]=%f\n",i,a+i*h,i,yy[i-1]);for(i=3;i<N;i++){x=a+i*h;yp=y+h*(55*f(x,y)-59*f(x-h,y2)+37*f(x-2*h,y1)-9*f(x-3*h,y0))/24;yc=y+h*(9*(*f)(x+h,yp)+19*(*f)(x,y)-5*(*f)(x-h,y2)+(*f)(x-2*h,y1))/24;printf("x[%d]=%f\ty[%d]=%f\n",i+1,x+h,i+1,yc);y0=y1;y1=y2;y2=y;y=yc;}}float f(float x,float y){return y*y;void Runge_Kutta(float(*f)(float x,float y),float a,float,float b,float y0,int N,float yy[]); void Adams (float a,float b,int N,float(*f)(float x,float y),float y0);float f(float x,float y);int main (void){float a=0,b=1,y0=1;int N=10;Adams(a,b,N,f,y0);}2、#include<stdio.h>void Runge_Kutta(float(*f)(float x,float y),float a,float b,float y0,int N,float yy[]){float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;for(i=1;i<=3;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/2,y+h*K3);y=y+h*(K1+2*K2+2*K2+2*K3+K4)/6;x=a+i*h;yy[i-1]=y;}}void Adams (float a,float b,int N,float(*f)(float x,float y),float y0)int i;float y1,y2,y,yp,yc,yy[3],h,x;printf("x[0]=%f\ty[0]=%f\n",a,y0);Runge_Kutta(f,a,b,y0,N,yy);y1=yy[0];y2=yy[1];y=yy[2];h=(b-a)/N;for(i=1;i<=3;i++)printf("x[%d]=%f\ty[%d]=%f\n",i,a+i*h,i,yy[i-1]);for(i=3;i<N;i++){x=a+i*h;yp=y+h*(55*f(x,y)-59*f(x-h,y2)+37*f(x-2*h,y1)-9*f(x-3*h,y0))/24;yc=y+h*(9*(*f)(x+h,yp)+19*(*f)(x,y)-5*(*f)(x-h,y2)+(*f)(x-2*h,y1))/24;printf("x[%d]=%f\ty[%d]=%f\n",i+1,x+h,i+1,yc);y0=y1;y1=y2;y2=y;y=yc;}}float f(float x,float y){return 0.1*(x*x*x+y*y);}void Runge_Kutta(float(*f)(float x,float y),float a,float,float b,float y0,int N,float yy[]); void Adams (float a,float b,int N,float(*f)(float x,float y),float y0);float f(float x,float y);int main (void){float a=0,b=1,y0=1;int N=10;Adams(a,b,N,f,y0);}。
计算方法上机答案
上海电力学院数值分析上机实验报告题目:数值分析上机实验报告学生姓名:11111111111学号:111111*********专业:11112013年12月30日数值计算方法上机实习题1. 设⎰+=105dx xx I nn , (1) 由递推公式n I I n n 151+-=-,从0I 的几个近似值出发,计算20I ; (2) 粗糙估计20I ,用nI I n n 51511+-=-,计算0I ;(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
(1) 解答:n=0,0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=⎰⎰⎰x d xdx x dx x x I nn这里可以用for 循环,while 循环,根据个人喜好与习惯:for 循环程序: While 循环程序: I=0.1823; I=0.1823; for n=1:20 i=1;I=(-5)*I+1/n; while i<21 End I=(-5)*I+1/i; I i=i+1; fprintf('I20=%f',I) end I = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009 (2) 粗略估计I 20: Mathcad 计算结果: for 循环程序: While 循环程序: >> I=0.007998; I=0.007998; >> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> II =0.0083(3) 算法误差分析:计算在递推过程中传递截断误差和舍入误差 第一种算法:(从1——>20)*000e I I=-***21111120115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===1x x 205x +⎛⎜⎜⎜⎠d 7.998103-⨯=误差放大了5n 倍,算法稳定性很不好; 第二种算法:(从20——>1)*n n ne I I =-***111111111()()555555n n n n n n nn e I I I I I I e n n ---=-=-+--+=-=0111...()55nne e e ===误差在逐步缩小,算法趋近稳定,收敛。
数值计算方法上机答案
0.23
22.322
b
54 240.236
29.304
117.818
2.991 1 13
3.907 0.017
3
1.007 2.1 0.5 61.705 0.91 4.5
1.006 6.3 1
12.17 4.918
[ 1.0, 2.0, 3.0, 4.5, 5.0, 21.8]
ans =
1.0000 -2.0000
2.9999 -4.0001
5.0000 -6.0008
17 1 4 3 1 2 3 7
2
10 1
7
2 1
1
4
1 1 8 2 5 2 1 1
3、
A
max1=0; for i=k:n
if abs(A(i,k))>max1 max1=abs(A(i,k));r=i;
end end if max1<1e-6
flag='failure';return; end if r>k
for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
[ -1.0, 34.211, -1.0, -2.1, 6.3, -1.7]
[
0, 0.5, 13.0, -0.5, 1.0, -1.5]
[ 4.501, 3.11, -3.907, -61.705, 12.17, 8.999]
[ 0.101, -0.812, -0.017, -0.91, 4.918, 0.1]
k x(k) 0 1.500000 1 0.784472 2 0.739519
数值计算方法上机实习题答案.doc
1.设I n 1 x ndx ,0 5 x( 1)由递推公式 I n 5I n 11,从 I 0的几个近似值出发,计算I 20;n解:易得: I 0 ln6-ln5=0.1823, 程序为:I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为: I 20= -3.0666e+010( 2)粗糙估计 I 20,用 I n 1 1I n 1 1 ,计算 I 0;5 5n0.0079 1 x 20 1 x 200.0095因为dx I 20dx 6 5所以取 I 20 1(0.0079 0.0095) 0.0087 2程序为: I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endII 0= 0.0083( 3)分析结果的可靠性及产生此现象的原因(重点分析原因 )。
首先分析两种递推式的误差;设第一递推式中开始时的误差为E0 I 0 I 0,递推过程的舍入误差不计。
并记 E n I n I n,则有 E n 5E n 1 ( 5) n E0。
因为 E20( 5) 20 E0 I 20,所此递推式不可靠。
而在第二种递推式中E0 1E1 (1)n E n,误差在缩小,5 5所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。
2.求方程e x10x 2 0 的近似根,要求x k 1x k 5 10 4,并比较计算量。
(1)在 [0, 1]上用二分法;程序: a=0;b=1.0;while abs(b-a)>5*1e-4c=(b+a)/2;if exp(c)+10*c-2>0b=c;else a=c;endendc结果: c =0.0903( 2)取初值x0 0,并用迭代 x k 1 2 e x ;10程序: x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果: x =0.0905(3)加速迭代的结果;程序: x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果: x =0.0995( 4)取初值x00 ,并用牛顿迭代法;程序: x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x;end x 结果: x =0.0905( 5) 分析绝对误差。
数值计算方法上机实习题NEW
数值计算方法上机实习题1. 设,(1) 由递推公式,从 , 出发,计算 ; 程序如下:function I=myhs(I0,n) if n>=1I=myhs(I0,n-1)*(-5)+1/n; elseif n==0 I=I0; end命令行窗口输入: I0=0.1822;I1=myhs(I0,20); I0=0.1823;I2=myhs(I0,20);运行结果:当I0=0.1822时,I20=-1.1593e+10。
当I0=0.1823时,I20= -2.0558e+9。
(2) ,20=10000I , 用,计算0I ; 程序如下:function I=myhs2(I20,n) if (n<20)I=myhs2(I20,n+1)*(-1)/5+1/(5*(n+1)); elseif n==20 I=I20; end命令行窗口输入:I20=0;I1=myhs2(I20,20); I20=10000;I2=myhs2(I20,20);运行结果:当I 20=0时,I 0=0.182321556793955。
当I 20=10000时,I 0= 0.182321556898812。
(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
根据上式,假设I n 的真值为I*,误差为e n ,即e=I*-I n 。
综合递推式,有e n =-5*e n-1,这意味着哪怕开始只有一点点误差,只要n 足够大,按照这种计算一步误差增长5倍的方式,所得的结果总是不可信的,因此整个算法是不稳定的。
根据上式,假设I n 的真值为I*,误差为e n ,即e=I*-I n 。
综合递推式,有e n-1. =(-1/5)*e n ,按照这种计算误差会以每步缩小到1/5的方式进行,根据(2)得到的结果而言,该算法是相对稳定的。
2. 求方程0210=-+x e x 的近似根,要求41105-+⨯<-k k x x ,并比较计算量。
数值计算方法上机实验题
数值计算方法上机实验实验内容:1.要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.2.给定积分dx e x⎰31和dx x ⎰311 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1)变步长梯形法; 2)变步长 Simpson 法; 3) Romberg 方法.算法描述:1、复合梯形法:⎰=tdt t a t V 0)()( ))()(2)((211∑-=++=n k k n b f x f a f hT输入 被积函数数据点t,a. 输出 积分值.n T复合Simpson 法:⎰=tdt t a t V 0)()( ))()(2)(4)((6101121∑∑---=++++=n k n k k k n b f x f x f a f hS输入 被积函数f(x),积分区间[a,b]和n 输出 复合Simpson 积分值n S步1 .);()(;a x b f a f S nab h n ⇐-⇐-⇐ 步2 对n k ,,2,1 =执行).(2;2);(4;2x f S S hx x x f S S h x x n n n n +⇐+⇐+⇐+⇐步3 n n S hS ⨯⇐6步4 输出n SRomberg 积分法:根据已知数据对其进行多项式拟合得出p(x);f(x)⇐p(x); 输入 被积函数f(x),积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111⇐===+⇐-⇐k R C S b f a f hT a b h 步2 反复执行步3→步9. 步3 .2;0h a x S +⇐⇐ 步4 反复执行步5→步6. 步5 ;);(h x x x f S S +⇐+⇐步6 若x ≥b,则退出本层循环. 步7 执行.6316364;1511516;3134;2212212212212C C R S S C T T S S h T T -⇐-⇐-⇐+⇐步8 执行.1;;;;;2;2121212112+⇐⇐⇐⇐⇐⇐-⇐k k R R C C S S T T hh R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ⇐ 步11 输出.2n R2、变步长梯形算法:功能 求积分⎰ba)(dx x f ,允许误差为ε。
数值计算方法》习题答案
《数值计算方法》课后题答案详解吉 林 大 学第一章 习 题 答 案1. 已知(1)2,(1)1,(2)1f f f −===,求()f x 的Lagrange 插值多项式。
解:由题意知:()01201212001020211012012202121,1,2;2,1,1()()(1)(2)()()6()()(1)(2)()()2()()(1)(1)()()3(1)(2)(1)(2)()2162nj j j x x x y y y x x x x x x l x x x x x x x x x x l x x x x x x x x x x l x x x x x x x x L x y l x ==−=====−−−−==−−−−+−==−−−−−+−==−−−−+−==×+×−∴∑()2(1)(1)131386x x x x +−+×=−+2. 取节点01210,1,,2x x x ===对x y e −=建立Lagrange 型二次插值函数,并估计差。
解11201201210,1,;1,,2x x x y y e y e −−======1)由题意知:则根据二次Lagrange插值公式得:02011201201021012202110.510.520.51()()()()()()()()()()()()()2(1)(0.5)2(0.5)4(1)(224)(43)1x x x x x x x x x x x x L x y y y x x x x x x x x x x x x x x x x e x x e e e x e e x −−−−−−−−−−−−=++−−−−−−=−−+−−−=+−+−−+22)Lagrange 根据余项定理,其误差为(3)2210122()1|()||()||(1)(0.5)|3!61max |(1)(0.5)|,(0,1)6()(1)(0.5),()330.5030.2113()61()0.2113(0.21131)(0.21130.5)0.008026x f R x x e x x x x x x t x x x x t x x x x t x R x ξξωξ−+≤≤==−−≤−−∈′=−−=−+=−==≤××−×−=∴取 并令 可知当时,有极大值3. 已知函数y =在4, 6.25,9x x x ===处的函数值,试通过一个二次插值函数求的近似值,并估计其误差。
数值分析上机题答案
数值分析上机题答案【篇一:数值分析上机试题对应参考答案】么是近似值x* 有效数字?若近似值x*的误差限是某一位的半个单位,该位到x*的第一位非零数字共有n位,就说x*有n位有效数字。
它可表示为2、数值计算应该避免采用不稳定的算法,防止有效数字的损失. 因此,在进行数值运算算法设计过程中主要注意什么?(1)简化计算过程,减少运算次数;(2)避免两个相近的数相减;(3)避免除数的绝对值远小于被除数的绝对值;(4)防止大数“吃掉”小数的现象;(5)使用数值稳定的算法,设法控制误差的传播。
3、写出“n 阶阵a 具有n 个不相等的特征值”的等价条件(至少写3 个)(1)|a|不为零(2)n阶矩阵a的列或行向量组线性无关(3)矩阵a为满秩矩阵(4)n阶矩阵a与n阶可逆矩阵b等价4、迭代法的基本思想是什么?就是用某种极限过程去逐步逼近线性方程组精确解得方法。
其基本思想为:先任取一组近似解初值x0,然后按照某种迭代原则,由x0计算新的近似解x1,以此类推,可计算出x2,x3,…xk,。
,如果{x}收敛,则取为原方程组的解。
5、病态线性方程组的主要判断方法有哪些?(1)系数矩阵的某两行(列)几乎近似相关(2)系数矩阵的行列式的值很小(3)用主元消去法解线性方程组时出现小主元(4)近似解x*已使残差向量r=b-ax*的范数很小,但该近似解仍不符合问题要求。
6、lagrange 插值的前提条件是什么?并写出二次lagrange 插值的基函数。
1,j?i?(x)? 前提条件是:l i ,j?0,1,2?,n.?ij0,j?i?二次lagrange 插值的基函数: (x?x)(x?x)12??lx0(xx)(xx) 0?10?2 (x?x)(x?x)02?? lx1(xx)(xx)1?01?2(x?x)(x?x)01?? lx2(x?x)(x?x)20217、什么是数值积分的代数精度?如果某一个求积公式对于次数不超过m的多项式均能准确地成立,但对于m+1次多项式就不准确成立,则称该求积公式具有m次代数精度(或代数精确度)。
数值计算方法上机实习题考证
---------------------------------------------------此文档包含我们计算方法的经典算法包含(数值计算方法上机实习题)1. 设⎰+=105dx x x I nn , (1) 由递推公式nI I n n 151+-=-,从0I 的几个近似值出发,计算20I ; (2) 粗糙估计20I ,用n I I n n 51511+-=-,计算0I ; (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
(1) 解答:n=0,0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=⎰⎰⎰x d x dx x dx x x I nn 这里可以用for 循环,while 循环,根据个人喜好与习惯:for 循环程序: While 循环程序:I=0.1823; I=0.1823;for n=1:20 i=1;I=(-5)*I+1/n; while i<21End I=(-5)*I+1/i;I i=i+1;fprintf('I20=%f',I) endI = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009(2) 粗略估计I 20: Mathcad 计算结果: for 循环程序: While 循环程序: >> I=0.007998; I=0.007998;>> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n);>> I n=n+1;I =0.0083 end>> II =0.0083(3) 算法误差分析:计算在递推过程中传递截断误差和舍入误差第一种算法:(从1——>20)01x x 205x +⎛⎜⎜⎜⎠d 7.998103-⨯=*000e I I =-***21111120115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===误差放大了5n 倍,算法稳定性很不好;第二种算法:(从20——>1)*n n ne I I =- ***111111111()()555555n n n n n n n n e I I I I I I e n n ---=-=-+--+=-=0111...()55n n e e e === 误差在逐步缩小,算法趋近稳定,收敛。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
1. 设⎰+=105dx xx I nn , (1) 由递推公式nI I n n 151+-=-,从0I 的几个近似值出发,计算20I ; 解:易得:0I =ln6-ln5=0.1823, 程序为:I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为:20I = -3.0666e+010 (2) 粗糙估计20I ,用nI I n n 515111+-=--,计算0I ; 因为 0095.056 0079.01020201020≈<<≈⎰⎰dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为:I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
首先分析两种递推式的误差;设第一递推式中开始时的误差为000I I E '-=,递推过程的舍入误差不计。
并记nn n I I E '-=,则有01)5(5E E E n n n -==-=- 。
因为=20E 20020)5(I E >>-,所此递推式不可靠。
而在第二种递推式中n nE E E )51(5110-==-= ,误差在缩小,所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。
2. 求方程0210=-+x e x的近似根,要求41105-+⨯<-k k x x ,并比较计算量。
(1) 在[0,1]上用二分法; 程序:a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果:c =0.0903(2) 取初值00=x ,并用迭代1021x k e x -=+;程序:x=0; a=1;while abs(x-a)>5*1e-4 a=x;x=(2-exp(x))/10; end x结果:x =0.0905(3) 加速迭代的结果; 程序:x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果:x =0.0995(4) 取初值00=x ,并用牛顿迭代法; 程序:x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x结果: x =0.0905(5) 分析绝对误差。
solve('exp(x)+10*x-2=0')3.钢水包使用次数多以后,钢包的容积增大,数据如下:试从中找出使用次数和容积之间的关系,计算均方差。
(注:增速减少,用何种模型) 设y=f(x)具有指数形式xbaey =(a>0,b<0)。
对此式两边取对数,得xba y 1ln ln +=。
记A=lna ,B=b ,程序:t=[0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625];z=[1.8594 2.1041 2.2597 2.2513 2.2721 2.3026 2.2956 2.3016 2.3504 2.3599 2.3609 2.3795 2.3609 2.3888 2.3758]; polyfit(t,z,1) 结果:ans = -1.1107 2.4578由此可得 A=2.4578,B=-1.1107,6791.11==Ae a ,b=B=-1.1107方程即为xey 1107.16791.11-=计算均方差编程:x=[2:16]; y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76]; f(x)=11.6791*exp( -1.1107./x); c=0;for i=1:15 a=y(i); b=x(i);c=c+(a-f(b))^2; endaverge=c/15 结果:averge =0.05944.设⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛----------------=410100141010014101101410010141001014A ,⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=625250b ,b x =A 分析下列迭代法的收敛性,并求42110-+≤-kk x x 的近似解及相应的迭代次数。
(1) JACOBI 迭代; 程序:function y=jacobi(a,b,x0) D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1); B=D\(L+U); f=D\b;y=B*x0+f;n=1;while norm (y-x0)>1e-4 x0=y;y=B*x0+f;n=n+1; end y n以文件名jacobi.m 保存。
程序: a=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]'; jacobi(a,b,x0);运行结果为:y =1.00002.00001.00002.00001.00002.0000n =28(2)GAUSS-SEIDEL迭代;程序:function y=seidel(a,b,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>10^(-4)x0=y;y=G*x0+f;n=n+1;endyn以文件名deisel.m保存。
程序:a=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]';x0=[0 0 0 0 0 0]';jacobi(a,b,x0);运行结果为:y =1.00002.00001.00002.00001.00002.0000n =15(3) SOR 迭代(95.0,95.1,334.1=ω)。
程序:function y=sor(a,b,w,x0) D=diag(diag(a)); U=-triu(a,1); L=-tril(a,-1);lw=(D-w*L)\((1-w)*D+w*U); f=(D-w*L)\b*w; y=lw*x0+f;n=1;while norm(y-x0)>10^(-4) x0=y;y=lw*x0+f;n=n+1; end y n以文件名sor.m 保存。
程序:a=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]'; x0=[0 0 0 0 0 0]'; c=[1.334 1.95 0.95]; for i=1:3 w=c(i); sor(a,b,w,x0); end运行结果分别为: y =1.00002.0000 1.0000 2.0000 1.0000 2.0000 n =13 y =1.00002.0000 1.0000 2.0000 1.0000 2.0000 n =241 y =1.00002.0000 1.0000 2.0000 1.0000 2.0000 n =175.用逆幂迭代法求⎪⎪⎪⎭⎫ ⎝⎛=111123136A 最接近于11的特征值和特征向量,准确到310-。
程序:function [mt,my]=maxtr(A,p,ep) n=length(A); B=A-p*eye(n); v0=ones(n,1); k=1; v=B*v0;while abs(norm(v,inf)-norm(v0,inf))>ep %norm(v-v0)>ep k=k+1;q=v;u=v/norm(v,inf) v=B*u; v0=q; endmt=1/norm(v,inf)+p my=u主界面中输入:A=[1 -2 -3]; maxtr(A,11,0.001) 结果为: 特征值: mt =11.0919特征向量: my =0.3845 -1.0000 0.73066.用经典R-K 方法求解初值问题(1)⎩⎨⎧-+-='++-='x x y y y x y y y sin 2cos 22sin 22212211,]10,0[∈x ,⎩⎨⎧==3)0(2)0(21y y ; 程序:function ydot=lorenzeq(x,y)ydot=[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)] 以文件民lorenzeq.m 保存。
主窗口输入:[x,y]=ode45('lorenzeq',[0:10],[2;3]) 运行结果为: x =0 1 2 3 4 5 6 7 8 9 10 y =2.00003.0000 1.5775 1.2758 1.1802 -0.1457 0.2406 -0.8903 -0.7202 -0.6170 -0.9454 0.2971 -0.2745 0.9652 0.6589 0.7557 0.9901 -0.1449 0.4124 -0.9109 -0.5440 -0.8389 (2)⎩⎨⎧-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212211,]10,0[∈x ,⎩⎨⎧==3)0(2)0(21y y 。
和精确解⎩⎨⎧+=+=--xe x y xe x y xx cos 2)(sin 2)(21比较,分析结论。
程序:function ydot=lorenzeq1(x,y)ydot=[-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)]; 以文件名lorenzeq1.m 保存。