A4 Scripts
Specifications of Functions that you must write...
function X = GenSylvester(S,R,Q,T,B) % S and Q are m-by-m upper triangular matrices % R and T are n-by-n upper triangular matrices % B is an m-by-n matrix % X solves the linear system SXR - QXT = B. % It is assumed that this system is nonsingular function Z = Evolve(v,k) % v is a column vector and k is a nonzero integer. % Define n = length(v)-1 and let P be the (n+1)-by-(n+1) % plasmid transition matrix. % Z is a (n+1)-by-|k| matrix % If k is positive then Z(:,i) = (P^i)v, i = 1:k. % If k is negative, then (P^i)Z(:,i) = v, i = 1:|k|. function [d,u,v] = InsideDist(x,y,theta) % Assume that x and y are column n vectors % and that theta is a scalar that satisfies 0 <= theta < 2pi. % % Let R be the ray that originates at (0,0) and makes an angle theta with the % positive x-axis. % % Let P be the polygon defined by connecting % (x(1),y(1)),...,(x(n),y(n)),(x(1),y(1)) % Assume that % (i) (0,0) is outside P, % (ii) only adjacent polygon edges intersect, % (iii) the ray does not go through any polygon vertex and is % not parallel to any edge. % % d is the length of that part of R that is inside P. % % If R intersects the polygon m times, then u and v are column m-vectors % that define the intersection points (u(1),v(1)),...,(u(m),v(m)). It is assumed % that these points are indexed in order of increasing distance from (0,0). % % If R does not intersect the polygon, then u and v are empty vectors.
A4A
% A4A
clc
randn('seed',0)
% Make sure things can be OK even if some of the triangular
% matrices are singular.
M1 = [1 2; 0 3]; M2 = [ 10 8; 0 5] ;
M3 = [ 0 1 ; 0 1]; M4 = [4 1 ; 0 13 ];
B = randn(2,2);
X = GenSylvester(M1,M2,M3,M4,B)
Residual = norm(M1*X*M2-M3*X*M4-B,1)
X = GenSylvester(M2,M3,M4,M1,B)
Residual = norm(M2*X*M3-M4*X*M1-B,1)
X = GenSylvester(M3,M4,M1,M2,B)
Residual = norm(M3*X*M4-M1*X*M2-B,1)
X = GenSylvester(M4,M1,M2,M3,B)
Residual = norm(M4*X*M1-M2*X*M3-B,1)
% A larger example...
disp(' ' )
m = 11;
n = 6;
S = triu(randn(m,m)); Q = triu(randn(m,m));
R = triu(randn(n,n)); T = triu(randn(n,n));
B = randn(m,n);
X = GenSylvester(S,R,Q,T,B)
Normalized_Residual = norm(S*X*R - Q*X*T - B,1)/norm(X,1)
% Benchmark
disp(' ' )
nRepeat = 100; % Repeat factor to address clock granularity
m = 100;
% Examine time relative to matrix-matrix multiplication
A = randn(m,m);
tic
for k=1:nRepeat
C = A*A;
end
MatMultTime = toc;
S = triu(randn(m,m))+10*eye(m,m); Q = triu(randn(m,m))-5*eye(m,m);
R = triu(randn(m,m))+10*eye(m,m); T = triu(randn(m,m))-5*eye(m,m);
B = randn(m,m);
tic
for k=1:nRepeat
X = GenSylvester(S,R,Q,T,B);
end
GenSylvesterTime = toc;
% NOTE: Optimize GenSylvester for speed.
ratio = GenSylvesterTime/MatMultTime;
disp(sprintf('GenSylvesterTime/MatMultTime = %6.3f (m = %3d)',ratio,m))
A4B
% A4B
close all
bigscreen
rand('seed',0)
% Let's start with a random distribution and see how it evolves...
n = 20;
nDivides = 50;
v = rand(n+1,1); v = v/sum(v);
Z = Evolve(v,nDivides);
shg
for k=0:nDivides
if k==0
bar((0:n)',v);
else
bar((0:n)',Z(:,k))
end
title(sprintf('Distribution of type A Plasmids after %1d subdivisions',k))
axis([-1 n+2 0 .5])
set(gca,'XTick',0:n)
pause(.2)
end
text(5,.45,'Final Numerical Distribution of type A Plasmids')
for j=1:n+1
text(10,.4-.015*j,sprintf('%8.6f',Z(j,nDivides)))
end
pause
% Now let's explore some accuracy issues...
n = 11;
nStepsBack = 8;
v = rand(n+1,1); v = v/sum(v);
% Think of v as the current iA distribution
Z = Evolve(v,-nStepsBack);
shg
for j=1:nStepsBack
w = Z(:,j)/sum(Z(:,j));
% w is the predicted iA distribution j "timesteps" back in time
% Let's march forward j timesteps from this point...
Z0 = Evolve(w,j);
v0 = Z0(:,j);
% v0 should equal v, ...
plot((0:n)',abs(v0-v))
axis([-1 n+2 -inf inf])
set(gca,'XTick',0:n)
title(sprintf('Recoverying the current i_{A} distribution. ( j = %3.0d)',j))
xlabel('i_{A}')
ylabel('Absolute Error')
pause
end
A4C
% A4C
x = [ 1.3034 5.4825 2.6200 4.5121 5.0211 ...
3.2414 5.4866 6.4727 7.6891 7.3329 ...
8.2069 6.0548 6.2206 4.1411 4.3822 ...
5.3315 3.4781 3.4781 1.6849 2.7096 ...
0.4041 2.3178 -0.6507 -0.5089 1.4048 ...
4.4637 3.9062 -0.5842 -1.4281 0.9860 ...
2.2065 2.1161 -0.8524 1.1216 -1.4400 ...
0.1211 -1.9885 -1.0271 -2.1874 -0.2074 ...
0.7269 0.7269 0.3050 0.4406 0.6365 ...
0.7871 1.2995 1.2995 0.9378 0.9981]';
y = [ 1.1322 -0.7623 1.5240 1.4486 2.4130 ...
3.5281 3.6034 2.6089 2.3829 4.5829 ...
6.0445 5.5925 6.9034 7.4610 6.1349 ...
5.3664 5.0349 5.8938 7.3705 5.4267 ...
5.6979 4.7938 4.8938 6.6925 7.8527 ...
8.0637 9.2842 9.0130 4.2062 3.7692 ...
4.1158 3.3322 3.3774 2.7445 2.4884 ...
1.9760 1.7048 1.3884 0.6952 1.1473 ...
2.4281 1.6897 0.6952 0.2432 -0.5390 ...
0.3185 -0.4130 0.3486 0.7103 1.4938]';
close all
for theta = linspace(0,pi,19)
% Draw the polygon and the ray.
fill(x,y,'m')
hold on
plot([0 20*cos(theta)],[0 20*sin(theta)],'--')
[d,u,v] = InsideDist(x,y,theta);
title(sprintf('Inside Distance = %6.3f theta = %4.0f (degrees)',d,180*theta/pi))
for k=1:length(u)/2
% This partof the ray is inside P...
plot(u(2*k-1:2*k),v(2*k-1:2*k),'k')
end
hold off
axis([-3 10 -1 10])
pause
end