Contents
Initialize a 2D matrix
some 2D matrix operations
some ODE solutions
conjugate gradient and some timing
Other links
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 % stiff1
Links 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