用Romberg方法求解积分.docx

1、用Romberg方法求解积分1.用Romberg方法求解积分，要求误差不超过解：Romberg.m文件：function I, step = Romberg(f, a, b,EPS)% Romberg.m 是用龙贝格公式求积分% f 为被积函数% EPS 为积分结果精度% a,b 为积分区间的上下限% I 为积分结果；step 为积分的子区间数m = 1k = 0Er = 0.1H =b-aS = zeros(1, 1)S(1, 1) = (H/2) * (subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)while Er E

2、PS 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)/(4n-1) end Er = abs(S(k+1, n+1)-S(k, n)endI = S(k+1, k+1)step = k命令：clearclcformat shorta = 0; b = 0.8;

3、 EPS = 1e-2;I, step = Romberg(x(1/2), a, b, EPS)计算结果：m = 1k = 0Er = 0.1000H = 0.8000S = 0S = 0.3578k = 1f1 = 0H = 0.4000x = 0.4000f1 = 0.6325S = 0.3578 0.4319m = 2S = 0.3578 0 0.4319 0.4566Er = 0.0988k = 2f1 = 0H = 0.2000x = 0.2000f1 = 0.4472x = 0.6000f1 = 1.2218S = 0.3578 0 0.4319 0.4566 0.4603 0m

4、= 4S = 0.3578 0 0.4319 0.4566 0.4603 0.4698S = 0.3578 0 0 0.4319 0.4566 0 0.4603 0.4698 0.4707Er = 0.0141k = 3f1 = 0H = 0.1000x = 0.1000f1 = 0.3162x = 0.3000f1 = 0.8640x = 0.5000f1 = 1.5711x = 0.7000f1 = 2.4077S = 0.3578 0 0 0.4319 0.4566 0 0.4603 0.4698 0.4707 0.4709 0 0m = 8S = 0.3578 0 0 0.4319 0

5、.4566 0 0.4603 0.4698 0.4707 0.4709 0.4745 0S = 0.3578 0 0 0.4319 0.4566 0 0.4603 0.4698 0.4707 0.4709 0.4745 0.4748S = 0.3578 0 0 0 0.4319 0.4566 0 0 0.4603 0.4698 0.4707 0 0.4709 0.4745 0.4748 0.4748Er = 0.0042I = 0.4748step = 3I = 0.4748step = 32.设方程组试用Jacobi迭代法求解此方程，当时终止迭代。解：Jacobi.m文件：function

6、Jacobi(A, b, max, eps) %max为最大迭代次数，eps为容许误差n = 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

7、 + 1 if sum(abs(x1 - x) = max fprintf(The Method is disconvergentn) break end x = x1endif k max for i = 1:n fprintf( x %d = %fn,i,x1(i) ) endend命令：clearclcformat shortA = 5 2 1; -1 4 2; 2 -3 10;b = -12 20 3;max = 100;eps = 1e-5Jacobi(A, b, max, eps)计算结果：i = 1A = 5 2 1 -1 4 2 2 -3 10b = -12 20 3D = 5

8、 0 0 0 4 0 0 0 10L = 0 0 0 1 0 0 -2 3 0U = 0 -2 -1 0 0 -2 0 0 0D0 = 0.2000 0 0 0 0.2500 0 0 0 0.1000x0 = 0 0 0B = 0 -0.4000 -0.2000 0.2500 0 -0.5000 -0.2000 0.3000 0f = -2.4000 5.0000 0.3000x = -2.4000 5.0000 0.3000x0 = -2.4000 5.0000 0.3000i = 2x = -4.4600 4.2500 2.2800x0 = -4.4600 4.2500 2.2800i =

9、 3x = -4.5560 2.7450 2.4670x0 = -4.5560 2.7450 2.4670i = 4x = -3.9914 2.6275 2.0347x0 = -3.9914 2.6275 2.0347i = 5x = -3.8579 2.9848 1.8865x0 = -3.8579 2.9848 1.8865i = 6x = -3.9712 3.0922 1.9670x0 = -3.9712 3.0922 1.9670i = 7x = -4.0303 3.0237 2.0219x0 = -4.0303 3.0237 2.0219i = 8x = -4.0139 2.9815

10、 2.0132x0 = -4.0139 2.9815 2.0132i = 9x = -3.9952 2.9900 1.9972x0 = -3.9952 2.9900 1.9972i = 10x = -3.9954 3.0026 1.9960x0 = -3.9954 3.0026 1.9960i = 11x = -4.0002 3.0031 1.9999x0 = -4.0002 3.0031 1.9999i = 12x = -4.0012 3.0000 2.0010x0 = -4.0012 3.0000 2.0010i = 13x = -4.0002 2.9992 2.0002x0 = -4.0

11、002 2.9992 2.0002i = 14x = -3.9997 2.9998 1.9998x0 = -3.9997 2.9998 1.9998i = 15x = -3.9999 3.0002 1.9999x0 = -3.9999 3.0002 1.9999i = 16x = -4.0000 3.0001 2.0000x0 = -4.0000 3.0001 2.0000i = 17x = -4.0000 3.0000 2.0000x0 = -4.0000 3.0000 2.0000i = 18x = -4.0000 3.0000 2.0000x0 = -4.0000 3.0000 2.0000i = 19x = -4.0000 3.0000 2.0000x0 = -4.0000 3.0000 2.0000i = 20x = -4.0000 3.0000 2.000x = -4.0000 3.0000 2.0000i = 20