计算方法上机实践一二四.docx

1、计算方法上机实践一二四计算方法上机课实验报告第一次实验一 实验目的：验证和比较各种递推公式的稳定性，从而加深对算法的数值稳定性的认识。二 实验课题：递推计算的稳定性 三 采用的具体算法实验一：%其中a=1, b=-（5108），c=5108。a=1;b=-(5*108+1);c=5*108;x1=(-b+sqrt (b2-4*a*c)/(2*a)x2=(-b-sqrt (b2-4*a*c)/(2*a)X1=-(b-sqrt(b2-4*a*c)/(2*a)X2=c/(a*X1)实验二：function try_stableglobal n aN = 20; % 计算 N 个值a =0.05; %

2、-% 方案I 用递推公式 % I (k) = - a*I (k-1) + 1/k I0 = log (a+1)/a); % 初值I = zeros (N,1); I (1) = -a*I0+1;for k = 2:NI (k) = - a*I (k-1) + 1/k;end %- %方案II 用递推公式 % I (k-1) = ( - I(k) + 1/k ) / a II = zeros (N,1);if a = N/(N+1) II (N) = 0.5*(1/(a+1)*(N+1)+1/(a*(N+1);else II (N) = 0.5*(1/(a+1)*(N+1)+1/N);end f

3、or k = N:-1:2 II (k-1) = ( - I(k) + 1/k ) / a; end%-% 调用 matlab 高精度数值积分命令 quadl 计算以便比较 III = zeros (N,1); for k = 1:N n = k; III (k) = quadl (f,0,1);end%-% 显示计算结果% clc fprintf (n 方案I结果 方案II结果 精确值) for k = 1:N, fprintf (n I (%2 .0f) %17 .7f %17 .7f %17 .7f,k,I (k),II (k),III (k) end% 把下面的这段保存为f.mfunc

4、tion y = f (x) % 定义函数global n a % 参量 n 为全局变量 y =x.n./(a+x); % 注意:这里一定要 点 运算 return实验结果：实验一：x1=500000000 x2=0X1=500000000 X2=1实验二：a=0.05 方案I结果 方案II结果 精确值 I ( 1) 0.8477739 0.8477739 0.8477739 I ( 2) 0.4576113 0.4576113 0.4576113 I ( 3) 0.3104528 0.3104528 0.3104528 I ( 4) 0.2344774 0.2344774 0.2344776

5、 I ( 5) 0.1882761 0.1882761 0.1882761 I ( 6) 0.1572529 0.1572529 0.1572529 I ( 7) 0.1349945 0.1349945 0.1349945 I ( 8) 0.1182503 0.1182503 0.1182503 I ( 9) 0.1051986 0.1051986 0.1051986 I (10) 0.0947401 0.0947401 0.0947401 I (11) 0.0861721 0.0861721 0.0861724 I (12) 0.0790247 0.0790247 0.0790247 I (

6、13) 0.0729718 0.0729718 0.0729718 I (14) 0.0677800 0.0677800 0.0677800 I (15) 0.0632777 0.0632777 0.0632777 I (16) 0.0593361 0.0593361 0.0593361 I (17) 0.0558567 0.0558567 0.0558567 I (18) 0.0527627 0.0527627 0.0527627 I (19) 0.0499934 0.0499934 0.0499934 I (20) 0.0475003 0.0476757 0.0475003-a=15 方案

7、I结果 方案II结果 精确值 I ( 1) 0.0319222 0.0319222 0.0319222 I ( 2) 0.0211673 0.0211673 0.0211673 I ( 3) 0.0158245 0.0158245 0.0158245 I ( 4) 0.0126326 0.0126326 0.0126326 I ( 5) 0.0105112 0.0105112 0.0105112 I ( 6) 0.0089993 0.0089993 0.0089993 I ( 7) 0.0078674 0.0078674 0.0078674 I ( 8) 0.0069883 0.0069883

8、 0.0069883 I ( 9) 0.0062862 0.0062862 0.0062859 I (10) 0.0057064 0.0057064 0.0057117 I (11) 0.0053136 0.0053136 0.0052337 I (12) 0.0036289 0.0036289 0.0048296 I (13) 0.0224896 0.0224896 0.0044838 I (14) -0.2659159 -0.2659159 0.0041831 I (15) 4.0554050 4.0554050 0.0039207 I (16) -60.7685756 -60.76857

9、56 0.0036893 I (17) 911.5874579 911.5874579 0.0034837 I (18) -13673.7563129 -13673.7563129 0.0032998 I (19) 205106.3973248 205106.3973248 0.0031344 I (20) -3076595.9098724 0.0030754 0.0029847-五 实验心得与体会：第二次实验实验目的：1 学习非线性方程 f (x) = 0 数值求根命令fzero；学习用割线法求近似根；学习绘制隐函数的图像2 学习非线性方程组 F (X) = 0 数值求根命令fsolve。实

10、验课题：求解非线性方程（组）的根采用的具体算法：实验一：function mysecant f = inline (x-exp (-x); a = 0.5; b = 0.6; delta = 0.001; epsilon = 0.00000000000001; max1 = 50; c,err,iter,yc = secant (f,a,b,delta,epsilon,max1)% % - function c,err,iter,yc = secant (f,a,b,delta,epsilon,max1)% % c,err,iter,yc = secant (f,a,b,delta,epsil

11、on,max1)% % 输入: f 连续函数% % a,b 迭代初值% % delta,epsilon 容差% % max1 最大迭代次数% % 输出: c 近似根% % err 误差% % iter 迭代次数% % yc = f (c)% for k = 1:max1 ya = feval (f,a); % ya = f (a) yb = feval (f,b); c = a-ya*(b-a)/(yb-ya); % 割线与 x 轴交点的横坐标 err = abs (c-b); % 相邻两次迭代的误差 relerr = err/(abs (c)+eps); % 相对误差,eps 是matlab

12、常数 (机器精度)约为1e-16 yc = feval (f,c); if (errdelta) | (relerrdelta) | (abs (yc)epsilon) % |是或 break end a = c; b = b; end iter = k;实验二：function implicit_function global p % 定义全局变量 n =101; x = linspace (-5,5,101); y = zeros (1,n); % 定义矩阵,初值是零,这是最常用的定义矩阵的方法 y0 = -4.6; % 第一次迭代初值 for k = 1:n p =x (k) ; y (

13、k) = fzero (fun,y0); y0 = y (k); end plot (x,y) % 作图 title (y3/(2+0.1*sin (x*y)+x2-4*x) % 加个标题% %- function z = fun (y) % 定义函数,这是最常用的定义函数的方式 global p x = p; z = y3/(2+0.1*sin (x*y)+x2-4*x;clf % 清图像x = -1.5:0.01:1.5; % 离散化,步长 0.01,这也是常用的方法y1 = -(x-3)/2; % 第一个方程求函数值,注意所有运算都是点运算,和常数相乘等就不用写点了y2 = (-2*x.

14、2+5).0.5; % 第二个方程求函数值,注意这里自变量与函数是颠倒的plot (x,y1,x,y2) % 作图,注意再把自变量与函数颠倒过来grid onaxis (-2 2 -1 4) % 限制横坐标与纵坐标的范围,这里要通过不断偿试来得到合适的范围% 通过作图发现在 (-1,2) 附近有一个根,调用 fsolve 求更精确的解clcX0 = -1,2; % 初值X = fsolve (myfun,X0)% - 定义函数组 myfun -function F = myfun (X)F = X (1)-3+2*X (2); 2*X (1)2+X (2)2-5; 实验结果：实验一：iter

15、= 5yc =7.7005e-013c = 0.5671err =0.0329iter = 6yc =-4.5519e-015实验二：实验：X = -0.8214 1.9107心得与体会：第三次实验实验目的：解线性方程组左除命令 的学习与矩阵计算有关的一些常用函数病态矩阵试验OR迭代法收敛速度受松驰因子的影响试验实验课题：具体算法：exp3_1% 【实验二】 x = Ab 与 x = inv (A)*b 在耗时方面的区别% 解方程组 Ax = b 时尽量不要使用 x = inv (A)*b,而要使用 x = Ab% 下面是二者耗时方面的区别实验 (把下面程序拷贝为新的M-文件,有些命令可不关心

16、其含义) rand (state,0); A = gallery (randsvd,200,2e13,2); % 产生条件数为2e13的200阶的随机矩阵 x = ones (200,1); % 设精确解为1,1,.,1 b = A*x; format long % 用15位小数显示% % % 求逆法 tic % 启动计时器 x1 = inv (A)*b; time1 = toc % 关闭计时器并显示耗时 error1 = norm (x-x1,inf) % 计算最大误差% % % 左除法 tic x2 = Ab; time2 = toc error2 = norm (x-x2,inf)实验结

17、果：time1 = 0.179745959332820error1 = 0.019363403320313time2 = 0.005831327724613error2 = 0.005499554430242n =13A 的条件数 = 365015202015855100.000000Gauss法求解结果 pcg法求解结果 1.00 1.00 1.00 1.00 1.01 1.01 0.81 1.00 2.72 1.00 -8.32 1.00 33.46 1.00 -74.16 1.00 117.92 1.00 -119.77 1.00 80.42 1.00 -29.09 1.00 6.00

18、0.99-n =14A 的条件数 = 4084634689031585800.000000Gauss法求解结果 pcg法求解结果 1.00 1.00 1.00 1.00 1.00 1.01 1.03 1.00 0.83 0.99 1.41 1.00 1.69 1.00 -6.86 1.00 26.48 1.00 -45.46 1.01 53.34 1.00 -35.21 1.00 15.14 1.00 -1.39 0.99-实验二取V=1:3:30c2 = 1.8063e+015心得与体会：第四次实验一 实验目的：1 学习一维插值命令；2 通过Lagrange 插值的Runge现象和Newton 插值的Runge现象演示增加对龙格现象的认识；3 学习三次样条的，通过三次样条插值的龙格现象的演示增加对样条插值优点的感性认识。 二 实验课题：龙格现象的发生、防止和插值效果的对比 实验结果：Lagrange插值样条插值