0001 function [Tinit,Tprop,E] = ngsetup_vpa(den,fs, ndigits)
0002
0003
0004
0005 digits(ndigits)
0006 fs = sym(fs);
0007 den = sym(den,'d');
0008 den=den';
0009 dt = vpa(1/fs);
0010
0011 n = length(den)-1;
0012
0013
0014
0015 m_a = vpa(zeros(n,n));
0016 for i = 1:n
0017 for j = 1:n
0018 if j == i+1
0019 m_a(i,j) = 1;
0020 end
0021 if i == n
0022 m_a(i,j) = -den(j);
0023 end
0024 end
0025 end
0026
0027
0028 a = m_a*dt;
0029 E = expm(a);
0030
0031
0032
0033 B = vpa(zeros(n,n));
0034 for i=1:n
0035 if rem(i,2) ~= 0
0036 j0 = (i+1)/2;
0037 s = (-1)^(j0+1);
0038 j = j0;
0039 for k=1:2:(n+1)
0040 d1 = den(k);
0041 d2 = s*d1;
0042 B(i,j) = d2;
0043 s = -s;
0044 j = j+1;
0045 end
0046 end
0047 if rem(i,2) == 0
0048 j0 = i/2+1;
0049 s = (-1)^j0;
0050 j = j0;
0051 for k=2:2:(n+1)
0052 d1 = den(k);
0053 d2 = s*d1;
0054 B(i,j) = d2;
0055 s = -s;
0056 j = j+1;
0057 end
0058 end
0059 end
0060
0061
0062 m_k = vpa(zeros(n,1));
0063 m_k(n) = 0.5;
0064 m_m = vpa(B\m_k);
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086 Cinit = vpa(zeros(n,n));
0087 for i=1:n
0088 for j=1:n
0089 if rem((i+j),2) == 0
0090 d1 = (-1)^((i-j)/2);
0091 d2 = (i+j)/2;
0092 Cinit(i,j) = d1 * m_m(d2);
0093
0094
0095 end
0096 end
0097 end
0098
0099
0100
0101 Tinit = ao.mchol(Cinit);
0102
0103
0104 N = n*(n+1)/2;
0105
0106 m_d = vpa(zeros(N));
0107 g = zeros(n);
0108 for i=1:n
0109 for j=1:n
0110 if i>=j
0111 g(i,j) = (i*i-i)/2+(j);
0112 else
0113 g(i,j) = (j*j-j)/2+(i);
0114 end
0115 end
0116 end
0117
0118 for i=1:n
0119 for j=i:n
0120 for k=1:n
0121 m_d(g(i,j),g(j,k)) = m_d(g(i,j),g(j,k)) + m_a(i,k);
0122 m_d(g(i,j),g(i,k)) = m_d(g(i,j),g(i,k)) + m_a(j,k);
0123 end
0124 end
0125 end
0126
0127
0128
0129 m_q = vpa(zeros(1,g(n,n)));
0130 for i=1:n
0131 for j=i:n
0132 if i==n
0133 m_q(g(i,j)) = (E(n,n))*(E(n,n))-1;
0134 else
0135 m_q(g(i,j)) = (E(i,n)*E(j,n));
0136 end
0137 end
0138 end
0139
0140 m_p = m_d\m_q';
0141
0142 Cprop = vpa(zeros(n));
0143 for i=1:n
0144 for j=1:n
0145 Cprop(i,j) = m_p(g(i,j));
0146 end
0147 end
0148
0149 Tprop = ao.mchol(Cprop);
0150 Tprop = double(Tprop);
0151 E = double(E);
0152 Tinit = double(Tinit);
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169 end
0170