用Romberg方法求解積分.doc

1. 用Romberg方法求解積分，要求誤差不超過解：Romberg.m文件：function I, step = Romberg(f, a, b,EPS)% Romberg.m 是用龍貝格公式求積分% f 為被積函數(shù)% EPS 為積分結(jié)果精度% a,b 為積分區(qū)間的上下限% I 為積分結(jié)果；step 為積分的子區(qū)間數(shù)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 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)/(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; EPS = 1e-2;I, step = Romberg(x(1/2), a, b, EPS)計(jì)算結(jié)果：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 = 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.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. 設(shè)方程組試用Jacobi迭代法求解此方程，當(dāng)時(shí)終止迭代。解：Jacobi.m文件：function Jacobi(A, b, max, eps) %max為最大迭代次數(shù)，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 + 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)計(jì)算結(jié)果：i = 1A = 5 2 1 -1 4 2 2 -3 10b = -12 20 3D = 5 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 = 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 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.0002 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 = 1