1、1. 用Romberg方法求解积分,要求误差不超过 解:Romberg.m文件: function [I, step] = Romberg(f, a, b,EPS) % Romberg.m 是用龙贝格公式求积分 % f 为被积函数 % EPS 为积分结果精度 % a,b 为积分区间的上下限 % I 为积分结果;step 为积分的子区间数 m = 1 k = 0 Er = 0.1 H =b-a S = zeros(1, 1) S(1, 1) = (H/2) * (subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym
2、f)),b)) while Er > EPS k = k + 1 f1 = 0 H = H/2 for i = 1:m x = a +H*(2*i-1) f1 = f1 + subs(sym(f),findsym(sym(f)),x) end S(k+1, 1) = S(k, 1)/2 + H*f1 m = 2 * m for n = 1:k S(k+1, n+1) = S(k+1, n) + (S(k+1, n)-S(k, n))/(4^n-1)
3、 end Er = abs(S(k+1, n+1)-S(k, n)) end I = S(k+1, k+1) step = k 命令: clear clc format short a = 0; b = 0.8; EPS = 1e-2; [I, step] = Romberg('x^(1/2)', a, b, EPS) 计算结果: m = 1 k = 0 Er = 0.1000 H = 0.8000 S = 0 S =
4、 0.3578 k = 1 f1 = 0 H = 0.4000 x = 0.4000 f1 = 0.6325 S = 0.3578 0.4319 m = 2 S = 0.3578 0 0.4319 0.4566 Er = 0.0988 k = 2 f1 = 0 H =
5、 0.2000 x = 0.2000 f1 = 0.4472 x = 0.6000 f1 = 1.2218 S = 0.3578 0 0.4319 0.4566 0.4603 0 m = 4 S = 0.3578 0 0.4319 0.4566 0.4603 0.4698 S = 0.35
6、78 0 0 0.4319 0.4566 0 0.4603 0.4698 0.4707 Er = 0.0141 k = 3 f1 = 0 H = 0.1000 x = 0.1000 f1 = 0.3162 x = 0.3000 f1 = 0.8640 x = 0.5000 f1 =
7、 1.5711 x = 0.7000 f1 = 2.4077 S = 0.3578 0 0 0.4319 0.4566 0 0.4603 0.4698 0.4707 0.4709 0 0 m = 8 S = 0.3578 0 0 0.4319 0.4566 0 0
8、4603 0.4698 0.4707 0.4709 0.4745 0 S = 0.3578 0 0 0.4319 0.4566 0 0.4603 0.4698 0.4707 0.4709 0.4745 0.4748 S = 0.3578 0 0 0 0.4319 0.4566 0 0 0.4
9、603 0.4698 0.4707 0 0.4709 0.4745 0.4748 0.4748 Er = 0.0042 I = 0.4748 step = 3 I = 0.4748 step = 3 2. 设方程组试用Jacobi迭代法求解此方程,,当时终止迭代。 解:Jacobi.m文件: function Jacobi(A, b, max, eps) %max为最大迭代次数,eps为容许误差 n
10、 length(A); x = zeros(n, 1); x1 = zeros(n, 1); k = 0; while 1 x1(1) = ( b(1) - A(1,2:n) * x(2:n,1) )/A(1,1) for i = 2:n-1 x1(i) = ( b(i) - A(i,1:i-1) * x(1:i-1,1) - A(i,i+1:n) * x(i+1:n,1))/A(i,i) end x1(n) = ( b(n) - A(n,1:n-1) * x(1:n-1,1) )/A(n,n) k = k + 1
11、 if sum(abs(x1 - x)) < eps fprintf('number = %d\n',k) break end if k >= max fprintf('The Method is disconvergent\n') break end x = x1 end if k < max for i = 1:n fprintf( 'x[ %d ] = %f\n',i,x1(i) ) end end 命令: clear c
12、lc format short A = [5 2 1; -1 4 2; 2 -3 10]; b = [-12 20 3]'; max = 100; eps = 1e-5 Jacobi(A, b, max, eps) 计算结果: i = 1 A = 5 2 1 -1 4 2 2 -3 10 b = -12 20 3 D = 5 0 0 0 4 0 0 0 10 L =
13、 0 0 0 1 0 0 -2 3 0 U = 0 -2 -1 0 0 -2 0 0 0 D0 = 0.2000 0 0 0 0.2500 0 0 0 0.1000 x0 = 0 0 0 B = 0 -0.4000 -0.2000 0
14、2500 0 -0.5000 -0.2000 0.3000 0 f = -2.4000 5.0000 0.3000 x = -2.4000 5.0000 0.3000 x0 = -2.4000 5.0000 0.3000 i = 2 x = -4.4600 4.2500 2.2800 x0 = -4.4600 4.2500 2.2800 i = 3 x =
15、4.5560 2.7450 2.4670 x0 = -4.5560 2.7450 2.4670 i = 4 x = -3.9914 2.6275 2.0347 x0 = -3.9914 2.6275 2.0347 i = 5 x = -3.8579 2.9848 1.8865 x0 = -3.8579 2.9848 1.8865 i = 6 x = -3.9712
16、 3.0922 1.9670 x0 = -3.9712 3.0922 1.9670 i = 7 x = -4.0303 3.0237 2.0219 x0 = -4.0303 3.0237 2.0219 i = 8 x = -4.0139 2.9815 2.0132 x0 = -4.0139 2.9815 2.0132 i = 9 x = -3.9952 2.9900
17、 1.9972 x0 = -3.9952 2.9900 1.9972 i = 10 x = -3.9954 3.0026 1.9960 x0 = -3.9954 3.0026 1.9960 i = 11 x = -4.0002 3.0031 1.9999 x0 = -4.0002 3.0031 1.9999 i = 12 x = -4.0012 3.0000 2.0010 x0
18、 -4.0012 3.0000 2.0010 i = 13 x = -4.0002 2.9992 2.0002 x0 = -4.0002 2.9992 2.0002 i = 14 x = -3.9997 2.9998 1.9998 x0 = -3.9997 2.9998 1.9998 i = 15 x = -3.9999 3.0002 1.9999 x0 = -3.999
19、9 3.0002 1.9999 i = 16 x = -4.0000 3.0001 2.0000 x0 = -4.0000 3.0001 2.0000 i = 17 x = -4.0000 3.0000 2.0000 x0 = -4.0000 3.0000 2.0000 i = 18 x = -4.0000 3.0000 2.0000 x0 = -4.0000 3.0000 2.0000 i = 19 x = -4.0000 3.0000 2.0000 x0 = -4.0000 3.0000 2.0000 i = 20 x = -4.0000 3.0000 2.000 x = -4.0000 3.0000 2.0000 i = 20






