Assignment 3 Scripts
P1
% P1 Examines Nearest1
% Plots the function that is closest to sin(2pi*x) on [0,1] that is
% orthogonal to 1,x,...,x^k for k=0:5.
% In this setting, two functions f and g are orthogonal on [0,1]
% if f(x(1))*g(x(1)) + ... + f(x(m))*g(x(m)) = 0 where x = linspace(0,1,200).
m = 200;
x = linspace(0,1,m)';
b = sin(2*pi*x);
plot(x,b)
s = 'Orthogonal to f(x) = 1';
A = ones(m,1);
for k=1:6
v = Nearest1(A,b);
subplot(2,3,k)
plot(x,b,'.',x,v)
title(s)
A = [A x.^k];
s = [s ', x^' int2str(k)];
end
P2
% P2 Examines LSSum
p = 10; m = 20; n = 5;
randn('seed',0)
A = cell(p,1), B = cell(p,1); c = cell(p,1); alfa = zeros(p,1);
for k = 1:p
A{k} = randn(m,n); B{k} = randn(m,n); c{k} = randn(m,1);
end
Soln = [];
for q = 1:8
alfa = [1;ones(p-1,1)/10^q];
[x,y] = LSSum(A,B,alfa,c);
Soln = [Soln [x;y]];
end
clc
disp(' q=1 q=2 q=3 q=4 q=5 q=6 q=7 q=8')
disp(Soln)
disp('As q increases block rows 2:p are "tuned out".')
P3
% P3 Examines MinPhi
rand('seed',0)
randn('seed',0)
tol = .001;
m = 150;
A = [cos(1:m)' sin(1:m)' linspace(0,5,m)'];
clc
format long
for tol = [.1 .01 .0001 .000001]
% Generate a ranndom translation vector and random 3-by-3 Q
v0 = randn(3,1);
theta = rand(3,1)*2*pi; c = cos(theta); s = sin(theta);
Q0 = [1 0 0; 0 c(1) s(1) ; 0 -s(1) c(1)]*[c(2) s(2) 0 ; -s(2) c(2) 0;0 0 1]*[1 0 0; 0 c(3) s(3); 0 -s(3) c(3)];
% Rotate A, translate, and add noise to get B
e = ones(m,1);
B = (A*Q0' - e*v0') + tol*randn(m,3);
% Recover A from B...
[Q,v] = MinPhi(A,B);
disp(sprintf('\n------------------------------------------------------------------------------\n\ntol = %10.3e',tol))
Q0 = Q0, Q = Q,
disp(sprintf('v0'' = [ %18.14f %18.14f %18.14f ]',v0'))
disp(sprintf('v'' = [ %18.14f %18.14f %18.14f ]',v'))
end
P4
% P4 Tests TrigLS
clc
[Sunrise Sunset] = GenData;
tau = (1:365)';
P = 365.25;
disp('Time-of-Sunrise Model...')
for n = 1:6
[alfa beta gamma Predict] = TrigLS(tau,Sunrise,P,n);
err = norm(Predict - Sunrise,inf);
disp(sprintf('\n\nn = %1d',n))
disp(sprintf('alfa = %10.6f',alfa))
disp(['beta = ' sprintf('%10.6f',beta)])
disp(['gamma = ' sprintf('%10.6f',gamma)])
disp(sprintf('|| Predict - Sunrise ||_inf = %8.5f',err))
end
disp(' ')
disp('------------------------------------------------------------------------')
disp(' ')
disp('Time-of-Sunset Model...')
for n = 1:6
[alfa beta gamma Predict] = TrigLS(tau,Sunset,P,n);
err = norm(Predict - Sunset,inf);
disp(sprintf('\n\nn = %1d',n))
disp(sprintf('alfa = %10.6f',alfa))
disp(['beta = ' sprintf('%10.6f',beta)])
disp(['gamma = ' sprintf('%10.6f',gamma)])
disp(sprintf('|| Predict - Sunset ||_inf = %8.5f',err))
end
% The following produces a table of Sunrise and Sunset times for Ithaca during 2003.
% Go to http://aa.usno.navy.mil/data/docs/RS_OneYear.html and enter 'Ithaca NY'
% and a similar table is produced.
disp(' ')
[alfa beta gamma Predict_Rise] = TrigLS(tau,Sunrise,P,6);
[alfa beta gamma Predict_Set] = TrigLS(tau,Sunset,P,6);
PrintData(Predict_Rise,Predict_Set)
function [Sunrise,Sunset] = GenData()
Times = [...
0735 1644 0720 1719 0642 1755 0549 1832 0502 1905 0432 1936 0432 1947 0458 1926 0531 1841 0603 1748 0639 1700 0715 1634;...
0736 1645 0719 1721 0640 1757 0547 1833 0500 1906 0431 1937 0433 1947 0459 1925 0532 1839 0604 1746 0640 1658 0716 1634;...
0736 1645 0718 1722 0639 1758 0546 1834 0459 1908 0431 1938 0434 1947 0500 1924 0533 1837 0605 1745 0642 1657 0717 1634;...
0736 1646 0717 1723 0637 1759 0544 1835 0458 1909 0430 1939 0434 1946 0501 1922 0534 1836 0606 1743 0643 1656 0718 1633;...
0736 1647 0716 1725 0635 1800 0542 1836 0456 1910 0430 1939 0435 1946 0502 1921 0535 1834 0607 1741 0644 1655 0720 1633;...
0735 1648 0715 1726 0634 1801 0540 1837 0455 1911 0430 1940 0435 1946 0503 1920 0536 1832 0608 1739 0645 1653 0720 1633;...
0735 1649 0713 1727 0632 1803 0539 1838 0454 1912 0429 1941 0436 1945 0504 1919 0537 1830 0609 1738 0647 1652 0721 1633;...
0735 1650 0712 1729 0630 1804 0537 1839 0453 1913 0429 1941 0437 1945 0505 1917 0538 1829 0610 1736 0648 1651 0722 1633;...
0735 1651 0711 1730 0629 1805 0535 1841 0451 1914 0429 1942 0437 1945 0506 1916 0539 1827 0612 1734 0649 1650 0723 1633;...
0735 1652 0710 1731 0627 1806 0534 1842 0450 1915 0429 1942 0438 1944 0507 1915 0540 1825 0613 1733 0650 1649 0724 1633;...
0734 1654 0708 1733 0625 1807 0532 1843 0449 1916 0429 1943 0439 1944 0508 1913 0541 1823 0614 1731 0652 1648 0725 1633;...
0734 1655 0707 1734 0624 1809 0530 1844 0448 1917 0428 1943 0440 1943 0509 1912 0542 1822 0615 1729 0653 1647 0726 1633;...
0734 1656 0706 1735 0622 1810 0529 1845 0447 1918 0428 1944 0440 1943 0510 1910 0543 1820 0616 1728 0654 1646 0727 1633;...
0733 1657 0704 1736 0620 1811 0527 1846 0446 1919 0428 1944 0441 1942 0512 1909 0544 1818 0617 1726 0655 1645 0727 1634;...
0733 1658 0703 1738 0619 1812 0526 1847 0445 1921 0428 1945 0442 1941 0513 1908 0545 1816 0618 1724 0657 1644 0728 1634;...
0732 1659 0702 1739 0617 1813 0524 1848 0444 1922 0428 1945 0443 1941 0514 1906 0547 1814 0620 1723 0658 1643 0729 1634;...
0732 1700 0700 1740 0615 1814 0522 1850 0443 1923 0428 1945 0444 1940 0515 1905 0548 1813 0621 1721 0659 1642 0730 1634;...
0731 1702 0659 1742 0613 1816 0521 1851 0442 1924 0428 1946 0445 1939 0516 1903 0549 1811 0622 1720 0700 1642 0730 1635;...
0731 1703 0657 1743 0612 1817 0519 1852 0441 1925 0428 1946 0445 1939 0517 1902 0550 1809 0623 1718 0701 1641 0731 1635;...
0730 1704 0656 1744 0610 1818 0518 1853 0440 1926 0429 1946 0446 1938 0518 1900 0551 1807 0624 1717 0703 1640 0731 1636;...
0730 1705 0654 1745 0608 1819 0516 1854 0439 1927 0429 1947 0447 1937 0519 1859 0552 1806 0626 1715 0704 1639 0732 1636;...
0729 1707 0653 1747 0606 1820 0515 1855 0438 1928 0429 1947 0448 1936 0520 1857 0553 1804 0627 1714 0705 1639 0732 1636;...
0728 1708 0651 1748 0605 1821 0513 1856 0437 1928 0429 1947 0449 1935 0521 1855 0554 1802 0628 1712 0706 1638 0733 1637;...
0727 1709 0650 1749 0603 1822 0512 1857 0437 1929 0430 1947 0450 1934 0522 1854 0555 1800 0629 1711 0707 1637 0733 1638;...
0727 1710 0648 1750 0601 1824 0510 1859 0436 1930 0430 1947 0451 1933 0523 1852 0556 1759 0630 1709 0709 1637 0734 1638;...
0726 1712 0647 1752 0559 1825 0509 1900 0435 1931 0430 1947 0452 1932 0524 1851 0557 1757 0632 1708 0710 1636 0734 1639;...
0725 1713 0645 1753 0558 1826 0507 1901 0435 1932 0431 1947 0453 1931 0525 1849 0558 1755 0633 1706 0711 1636 0734 1640;...
0724 1714 0644 1754 0556 1827 0506 1902 0434 1933 0431 1947 0454 1930 0526 1847 0559 1753 0634 1705 0712 1635 0735 1640;...
0723 1716 -1 -1 0554 1828 0504 1903 0433 1934 0432 1947 0455 1929 0527 1846 0600 1752 0635 1704 0713 1635 0735 1641;...
0722 1717 -1 -1 0553 1829 0503 1904 0433 1935 0432 1947 0456 1928 0529 1844 0602 1750 0637 1702 0714 1635 0735 1642;...
0721 1718 -1 -1 0551 1830 -1 -1 0432 1935 -1 -1 0457 1927 0530 1842 -1 -1 0638 1701 -1 -1 0735 1643];
Sunrise = reshape(Times(:,1:2:24),31*12,1);
Sunrise = Sunrise(find(Sunrise>0));
SunriseHr = floor(Sunrise/100);
Sunrise = SunriseHr + (Sunrise - 100*SunriseHr)/60;
Sunset = reshape(Times(:,2:2:24),31*12,1);
Sunset = Sunset(find(Sunset>0));
SunsetHr = floor(Sunset/100);
Sunset = SunsetHr + (Sunset - 100*SunsetHr)/60;
function PrintData(Sunrise,Sunset)
h = floor(Sunrise);
m = floor((Sunrise - h)*60);
Rises = [int2str(zeros(length(Sunrise),1)) int2str(100*h+m)];
Rises = [Rises(1:59,:); ' ';' ';' ';Rises(60:120,:);' ';...
Rises(121:181,:);' ';Rises(182:273,:);' '; Rises(274:334,:);' ';Rises(335:365,:)];
h = floor(Sunset);
m = floor((Sunset - h)*60);
Sets = [int2str(100*h+m)];
Sets = [Sets(1:59,:); ' ';' ';' ';Sets(60:120,:);' ';...
Sets(121:181,:);' ';Sets(182:273,:);' '; Sets(274:334,:);' ';Sets(335:365,:)];
M = [];
for k=1:31
if k<=9
M = [M ; '0' int2str(k) ' '];
else
M = [M ; int2str(k) ' '];
end
end
b = [ ' ';' '; ' '; ' ';' ';' ';' '];b = [b;b;b;b;' ';' ';' '];
disp(' ')
disp(' Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec')
disp('Day Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set Rise Set')
disp(' h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m h m')
for k=1:12
M = [M b b b Rises(1+(k-1)*31:k*31,:) b b Sets(1+(k-1)*31:k*31,:)];
end
disp(M)