Links to .m source file and .out output file, click to read test_mtx.m source file test_mtx_m.out output file
Links to .m source file and .png output, shown below stiff1.m source file using ode23 stiff1_m.png plot output % stiff1.m using ode23 % y1' = y2 y2' = (1-y1^2)*y2-y1 function stiff1() [t,y]=ode23(@vdp, [0 20], [2; 0]); plot(t,y(:,1),'-o',t,y(:,2),'-o') title('Solution of van der Pol Equation (\mu = 1) with ODE23'); xlabel('Time t'); ylabel('Solution y'); legend('y_1','y_2') function dydt = vdp(t, y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)]; end % vdp end % stiff1Links to .m source file and .png output, shown below stiff.m source file using ode15s stiff_m.png plot output % stiff.m using ode15s % y1' = y2 y2' = 1000*(1-y1^2)*y2-y1 function stiff() [t,y]=ode15s(@vdp1000, [0 3000], [2 0]); plot(t,y(:,1),'-o') title('Solution with(\mu = 1000) with ODE15'); xlabel('Time t'); ylabel('Solution y'); legend('y_1','y_2') function dydt = vdp1000(t, y) dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)]; end % vdp1000 end % stiff
Links to .m source file and .png output, shown below difeq1.m ODE definition, first order run_difeq1.m source file using ode45 difeq1_m.png plot output % difeq1.m simple ODE function dydt=difeq1(t,y) dydt = t*t + 2.0*t + 3.0; % run_difeq1.m simple first order t_range = [0 1]; [t,y] = ode45(@difeq1, t_range, [0]); size(t); size(y); for i=1:20 j=(i-1.0)/19.0; t2(i) = j; y2(i) = j*j*j/3.0 + j*j + 3.0*j; end plot(t, y, t2, y2, '--'); xlabel('time'); ylabel('y'); title('ODE solution and exact solution');
Links to .m source file and .png output, shown below difeq1.m ODE definition, first order run_difeq1.m source file using ode45 difeq1_m.png plot output % difeq2.m second order function dydt=difeq2(t,y) dydt = [y(2);-y(1)]; % dy^2 = -y y=sin(t) % run_difeq2.m second order t_range = [0 6.28]; [t,y] = ode45(@difeq2, t_range, [0 1.0]); size(t); size(y); for i=1:64 j=(i-1.0)/10.0; t2(i) = j; y2(i) = sin(j); end plot(t, y(:,1), t, y(:,2), '.', t2, y2, '--'); xlabel('time'); ylabel('y'); title('ODE solution and first, and exact solution');
![]()
Links to .m source file and .png output, shown below driver_cg.m source file to be timed test_driver_cg.m source file various n test_driver_cg_m.out timing output various n driver_cg4.png one plot output function driver_cg(N,fid) h = 1 / (N+1); x = [h : h : 1-h]; y = x; [X, Y] = ndgrid(x,y); F = (-2*pi^2) * (cos(2*pi*X).*(sin(pi*Y)).^2 + (sin(pi*X)).^2.*cos(2*pi*Y)); b = h^2 * F(:); clear X Y F; tic; A = setupA(N); tol = 1.0e-6; maxit = 99999; u = zeros(N^2,1); [u,flag,relres,iter,resvec] = pcg(A,b,tol,maxit,[],[],u); Uint = reshape(u, [N N]); % N.B.: Uint has only solutions on interior points timesec = toc; % append boundary to x, y, and to U: x = [0 x 1]; y = [0 y 1]; [X, Y] = ndgrid(x,y); U = zeros(size(X)); U(2:end-1,2:end-1) = Uint; % plot numerical solution: figure; H = mesh(X,Y,U); xlabel('x'); ylabel('y'); zlabel('u'); % compute and plot numerical error: Utrue = (sin(pi*X)).^2 .* (sin(pi*Y)).^2; E = U - Utrue; figure; H = mesh(X,Y,E); xlabel('x'); ylabel('y'); zlabel('u-u_h'); % compute L^inf norm of error and print: enorminf = max(abs(E(:))); fprintf(fid,'N = %5d\n', N); fprintf(fid,'tol = %10.1e, maxit = %d\n', tol, maxit); fprintf(fid,'flag = %1d, iter = %d, relres = %24.16e\n', flag, iter, relres); fprintf(fid,'h = %24.16e\n', h); fprintf(fid,'h^2 = %24.16e\n', h^2); fprintf(fid,'enorminf = %24.16e\n', enorminf); fprintf(fid,'C = enorminf / h^2 = %24.16e\n', (enorminf/h^2)); fprintf(fid,'wall clock time = %10.2f seconds\n', timesec); function A = setupA(N) I = speye(N); s = [-1*ones(1,N-1) 2*ones(1,N) -1*ones(1,N-1)]'; i = [2:N 1:N 1:N-1]'; j = [1:N-1 1:N 2:N ]'; T = sparse(i,j,s); A = kron(I,T) + kron(T,I); % test_driver_cg.m Gobberts driver_cg.m function test_driver_cg fid=fopen('test_driver_cg_m.out','w'); N=4; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=8; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=16; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=32; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=64; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=128; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=256; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=512; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); N=1024; fprintf(fid,'test_driver_cg running, N=%d \n', N); driver_cg(N,fid); fprintf(fid,'driver_cg ran, N=%d \n\n', N); fprintf(fid,'test_driver_cg finished \n'); end % test_driver_cg test_driver_cg running, N=4 N = 4 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 3, relres = 4.0645864755649299e-16 h = 2.0000000000000001e-01 h^2 = 4.0000000000000008e-02 enorminf = 1.1672672560746400e-01 C = enorminf / h^2 = 2.9181681401865993e+00 wall clock time = 0.18 seconds driver_cg ran, N=4 test_driver_cg running, N=8 N = 8 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 10, relres = 1.3833872844782684e-15 h = 1.1111111111111110e-01 h^2 = 1.2345679012345678e-02 enorminf = 3.9152467804694280e-02 C = enorminf / h^2 = 3.1713498921802370e+00 wall clock time = 0.02 seconds driver_cg ran, N=8 test_driver_cg running, N=16 N = 16 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 24, relres = 6.6499448293299502e-07 h = 5.8823529411764705e-02 h^2 = 3.4602076124567475e-03 enorminf = 1.1267467797972608e-02 C = enorminf / h^2 = 3.2562981936140836e+00 wall clock time = 0.02 seconds driver_cg ran, N=16 test_driver_cg running, N=32 N = 32 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 48, relres = 5.5636640914657905e-07 h = 3.0303030303030304e-02 h^2 = 9.1827364554637292e-04 enorminf = 3.0127765546470453e-03 C = enorminf / h^2 = 3.2809136680106321e+00 wall clock time = 0.02 seconds driver_cg ran, N=32 test_driver_cg running, N=64 N = 64 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 96, relres = 7.0188563830089490e-07 h = 1.5384615384615385e-02 h^2 = 2.3668639053254440e-04 enorminf = 7.7810619700791062e-04 C = enorminf / h^2 = 3.2874986823584220e+00 wall clock time = 0.04 seconds driver_cg ran, N=64 test_driver_cg running, N=128 N = 128 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 192, relres = 9.3340277731590903e-07 h = 7.7519379844961239e-03 h^2 = 6.0092542515473829e-05 enorminf = 1.9764834768531969e-04 C = enorminf / h^2 = 3.2890661538314050e+00 wall clock time = 0.24 seconds driver_cg ran, N=128 test_driver_cg running, N=256 N = 256 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 387, relres = 8.9244309654496858e-07 h = 3.8910505836575876e-03 h^2 = 1.5140274644582053e-05 enorminf = 4.9797438427923169e-05 C = enorminf / h^2 = 3.2890710107258974e+00 wall clock time = 1.94 seconds driver_cg ran, N=256 test_driver_cg running, N=512 N = 512 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 783, relres = 9.0704526218887027e-07 h = 1.9493177387914229e-03 h^2 = 3.7998396467669058e-06 enorminf = 1.2494278882391185e-05 C = enorminf / h^2 = 3.2881068792000065e+00 wall clock time = 14.70 seconds driver_cg ran, N=512 test_driver_cg running, N=1024 N = 1024 tol = 1.0e-06, maxit = 99999 flag = 0, iter = 1581, relres = 9.3989188600245957e-07 h = 9.7560975609756097e-04 h^2 = 9.5181439619274241e-07 enorminf = 3.1266064945967287e-06 C = enorminf / h^2 = 3.2848909483856881e+00 wall clock time = 142.16 seconds driver_cg ran, N=1024 test_driver_cg finished![]()
some timing, using ge rather than cg
Links to .m source file and .png output, shown below driver_ge.m source file to be timed test_driver_ge.m source file various n test_driver_ge_m.out timing output various n driver_ge4.png one plot output function driver_ge(N,fid) h = 1 / (N+1); x = [h : h : 1-h]; y = x; [X, Y] = ndgrid(x,y); F = (-2*pi^2) * (cos(2*pi*X).*(sin(pi*Y)).^2 + (sin(pi*X)).^2.*cos(2*pi*Y)); b = h^2 * F(:); clear X Y F; tic; A = setupA(N); u = A \ b; Uint = reshape(u, [N N]); % N.B.: Uint has only solutions on interior points timesec = toc; % append boundary to x, y, and to U: x = [0 x 1]; y = [0 y 1]; [X, Y] = ndgrid(x,y); U = zeros(size(X)); U(2:end-1,2:end-1) = Uint; % plot numerical solution: figure; plot3(X,Y,U); view(27.5,20); xlabel('x'); ylabel('y'); zlabel('u'); % compute and plot numerical error: Utrue = (sin(pi*X)).^2 .* (sin(pi*Y)).^2; E = U - Utrue; figure; plot3(X,Y,E); view(8.5,2); xlabel('x'); ylabel('y'); zlabel('u-u_h'); % compute L^inf norm of error and print: enorminf = max(abs(E(:))); fprintf(fid,'N = %5d\n', N); fprintf(fid,'h = %24.16e\n', h); fprintf(fid,'h^2 = %24.16e\n', h^2); fprintf(fid,'enorminf = %24.16e\n', enorminf); fprintf(fid,'C = enorminf / h^2 = %24.16e\n', (enorminf/h^2)); fprintf(fid,'wall clock time = %10.2f seconds\n', timesec); function A = setupA(N) I = speye(N); s = [-1*ones(1,N-1) 2*ones(1,N) -1*ones(1,N-1)]'; i = [2:N 1:N 1:N-1]'; j = [1:N-1 1:N 2:N ]'; T = sparse(i,j,s); A = kron(I,T) + kron(T,I); % test_driver_ge.m Gobberts driver_ge.m function test_driver_ge fid=fopen('test_driver_ge_m.out','w'); N=4; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=8; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=16; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=32; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=64; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=128; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=256; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=512; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); N=1024; fprintf(fid,'test_driver_ge running, N=%d \n', N); driver_ge(N,fid); fprintf(fid,'driver_ge ran, N=%d \n\n', N); fprintf(fid,'test_driver_ge finished \n'); end % test_driver_ge test_driver_ge running, N=4 N = 4 h = 2.0000000000000001e-01 h^2 = 4.0000000000000008e-02 enorminf = 1.1672672560746400e-01 C = enorminf / h^2 = 2.9181681401865993e+00 wall clock time = 0.69 seconds driver_ge ran, N=4 test_driver_ge running, N=8 N = 8 h = 1.1111111111111110e-01 h^2 = 1.2345679012345678e-02 enorminf = 3.9152467804693947e-02 C = enorminf / h^2 = 3.1713498921802099e+00 wall clock time = 0.02 seconds driver_ge ran, N=8 test_driver_ge running, N=16 N = 16 h = 5.8823529411764705e-02 h^2 = 3.4602076124567475e-03 enorminf = 1.1267474206801520e-02 C = enorminf / h^2 = 3.2563000457656393e+00 wall clock time = 0.02 seconds driver_ge ran, N=16 test_driver_ge running, N=32 N = 32 h = 3.0303030303030304e-02 h^2 = 9.1827364554637292e-04 enorminf = 3.0127943178556160e-03 C = enorminf / h^2 = 3.2809330121447653e+00 wall clock time = 0.01 seconds driver_ge ran, N=32 test_driver_ge running, N=64 N = 64 h = 1.5384615384615385e-02 h^2 = 2.3668639053254440e-04 enorminf = 7.7812147640843321e-04 C = enorminf / h^2 = 3.2875632378256299e+00 wall clock time = 0.03 seconds driver_ge ran, N=64 test_driver_ge running, N=128 N = 128 h = 7.7519379844961239e-03 h^2 = 6.0092542515473829e-05 enorminf = 1.9766136774113097e-04 C = enorminf / h^2 = 3.2892828205801603e+00 wall clock time = 0.33 seconds driver_ge ran, N=128 test_driver_ge running, N=256 N = 256 h = 3.8910505836575876e-03 h^2 = 1.5140274644582053e-05 enorminf = 4.9807274579038996e-05 C = enorminf / h^2 = 3.2897206786709465e+00 wall clock time = 0.80 seconds driver_ge ran, N=256 test_driver_ge running, N=512 N = 512 h = 1.9493177387914229e-03 h^2 = 3.7998396467669058e-06 enorminf = 1.2500832120454497e-05 C = enorminf / h^2 = 3.2898314883078901e+00 wall clock time = 3.54 seconds driver_ge ran, N=512 test_driver_ge running, N=1024 N = 1024 h = 9.7560975609756097e-04 h^2 = 9.5181439619274241e-07 enorminf = 3.1313407333755094e-06 C = enorminf / h^2 = 3.2898648580026446e+00 wall clock time = 15.39 seconds driver_ge ran, N=1024 test_driver_ge finished![]()
A simple test of conjugate gradient
Links to .m source file and .png output, shown below (one of several versions of conjugate gradient) conjgrad.m source file to be tested test_conjgrad.m test source test_conjgrad_m.out test output % conjgrad.m solve Ax=b for x by conjugate gradient method (initial x 0 OK) % A must be symmetric positive definite for convergence % added file output of debug function [x] = conjgrad(A, b, x, tol, maxit) debug = 2; r = b - A * x; p = r; rsold = r' * r; for i = 1:maxit Ap = A * p; alpha = rsold / (p' * Ap); x = x + alpha * p; r = r - alpha * Ap; rsnew = r' * r; if(debug > 1) disp('conjgrad.m iter, rsnew =') disp(i) disp(rsnew) end % if if sqrt(rsnew) < tol break; end % if p = r + (rsnew / rsold) * p; rsold = rsnew; end % end for if(debug > 0) disp('conjgrad final iter, rsnew =') disp(i) disp(rsnew) end % if end % conjgrad % test_conjgrad.m simple case 4 by 4 matrix % x computed using b = 1,2,3,4 function test_conjgrad() diary('test_conjgrad_m.out'); disp('test_conjgrad_m running with simeq and conjgrad') disp('solve Ax=b for x given A and b') A = [ 1 2 5 9 ; 2 6 9 10 ; 5 9 11 13 ; 9 10 13 19 ] b = [1 ; 2 ; 3 ; 4 ] % must be column x = A \ b; % solve A x = b disp('simeq solution x=') disp(x) disp('initial guess for conjgrad') x0 = [0 ; 0 ; 0 ; 0 ] xc = conjgrad(A, b, x0, 1.0e-10, 9); disp('conjgrad solution xc=') disp(xc) disp('check A det, eig') Ad = det(A) [v e] = eig(A); disp(e) disp('l u lu(A)') [l u] = lu(A) disp('ichol(A)') C = ichol(A) % fails because A is not sparse disp('test_conjgrad.m finished') end % test_conjgrad test_conjgrad_m running with simeq and conjgrad solve Ax=b for x given A and b A = 1 2 5 9 2 6 9 10 5 9 11 13 9 10 13 19 b = 1 2 3 4 simeq solution x= 0.1567 0.0896 0.1194 0.0075 initial guess for conjgrad x0 = 0 0 0 0 conjgrad.m iter, rsnew = 1 0.3117 conjgrad.m iter, rsnew = 2 0.0467 conjgrad.m iter, rsnew = 3 5.0755e-04 conjgrad.m iter, rsnew = 4 7.2548e-24 conjgrad final iter, rsnew = 4 7.2548e-24 conjgrad solution xc= 0.1567 0.0896 0.1194 0.0075 check A det, eig Ad = 134.0000 -3.5791 0 0 0 0 -0.3051 0 0 0 0 3.2617 0 0 0 0 37.6226 l u lu(A) l = 0.1111 0.2353 1.0000 0 0.2222 1.0000 0 0 0.5556 0.9118 -0.8472 1.0000 1.0000 0 0 0 u = 9.0000 10.0000 13.0000 19.0000 0 3.7778 6.1111 5.7778 0 0 2.1176 5.5294 0 0 0 1.8611 ichol(A) First argument must be sparse. C = ichol(A)
Last updated 2/20/2017