sparse/RescalingRandom.m at master · egoolish/sparse · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
interval = 10000;
Nrange = 125:125:1000;%:25:125;                    %What values of N we are testing
prange = 2;                       %What values of p we are testing
m=1;                                %Counter Variable
n=1;
for p=prange
    %figure(m)
    y_full = zeros(0, max(Nrange));%Nrange(end));
    for each=1:max(Nrange)%Nrange(end)
        y_full(each) = 2^each;
%         r = rand;
%         if(r<(1/3))
%             y_full(each)=(each-1)*p-1;
%         elseif(r>=(1/3) && (r<2/3))
%             y_full(each)=(each-1)*p+1;
%         else
%             y_full(each)=(each-1)*p;
%         end
    end
    for N=Nrange
        disp(N);
        y = y_full(1:N);
        range = 2^N.*pi;
%         if (N == 500) t = linspace(250,range,interval);
%         else t = linspace(250,392699100,interval); end
        t = linspace(0,range,interval);
        results = zeros(size(t));
        results2 = zeros(size(t));
        count = 1;
        for each=t
           disp(count);
           results(count) = linearqn(N, each, y);
           count = count+1;
        end
        for m = n:length(Nrange)
            figure(m);
            plot(t, results, 'DisplayName', ['qn(s) N=', num2str(N)]);
            xlim([125 5000])
            hold on;
            xlabel(['concurrent quad qn']);
            legend('show');
        end
%         figure(m);
%         plot(t, results, 'DisplayName', ['qn(s) N=', num2str(N)]);
%         %plot(t, results, 'DisplayName', 'qn(s)');
%         hold on;
%         %plot(t, results2, 'DisplayName', '(sin(ps)/ps)^2');
%         xlabel(['concurrent cubic qn']);
%         %xlabel(['2totheN qn vs sin(ps)/ps for N = ', num2str(N)]);
%         legend('show');
         %print(figure(m), '-dpdf', ['2totheNCompareforN=', num2str(N), '.pdf']);
         %m = m+1;
         %figure(m);
%          plot(t, results2 - results);
%          xlabel(['2totheN sin(ps)/ps - qn(s) for N = ', num2str(N)]);
%          print(figure(m), '-dpdf', ['2totheNDifferenceforN=', num2str(N), '.pdf']);
          %m = m+1;
          n=n+1;
    end
%     h = figure(m);
%     set(h,'PaperOrientation','landscape');
%     set(h,'PaperUnits','normalized');
%     set(h,'PaperPosition', [0 0 1 1]);
%     %print(h, '-dpdf', ['CubicConcurrent', num2str(Nrange), '.pdf']);
%     m = m+1;
    %f = fit(t, results, 'smoothingspline');
    %disp(f);
end
for m = 1:8
    h = figure(m);
     set(h,'PaperOrientation','landscape');
     set(h,'PaperUnits','normalized');
     set(h,'PaperPosition', [0 0 1 1]);
     print(h, '-dpdf', ['QuadTail', num2str(m), '.pdf']);
end
% close all;

function result=Qn(y, N, t)
    result = (1./(2.*pi)).*(1./(N+1)).*func_sum(y, t);
end

function result=linearqn(N, s, y)
    %result = ((4*pi)/(N^3+N^2+2)).*(Qn(y, N, (s/(N^3)))); %Cubic
    %result = ((6*pi)/(2*N*N+N+3)).*(Qn(y, N, (s/(N^2)))); %Quad
    result = (((2*pi)*(N+1))/(2^(N+2)+N-1)).*(Qn(y, N, (s*(N+1))/(2^(N+2)))); %2^N

end

function result = func_sum(y, t)
    result = 3;
    for each=y
        result = result + (sin((each+.5).*t))./(sin(.5.*t));
    end
end