Commit 8ded2d23 authored by Simon's avatar Simon
Browse files
parents ce2782a1 97eb3632
...@@ -72,6 +72,24 @@ else ...@@ -72,6 +72,24 @@ else
LDisc = L; LDisc = L;
end end
end end
if isa(A0,'function_handle')
A0Disc = eval_pointwise(A0,zdisc);
else
if size(A0,1)~= ndisc % constant matrix
A0Disc(1:ndisc,:,:) = repmat(shiftdim(A0,-1),ndisc,1,1);
else
A0Disc = A0;
end
end
if isa(H,'function_handle')
HDisc = eval_pointwise(H,zdisc);
else
if size(H,1)~= ndisc % constant matrix
HDisc(1:ndisc,:,:) = repmat(shiftdim(H,-1),ndisc,1,1);
else
HDisc = H;
end
end
n = size(LDisc,2); n = size(LDisc,2);
T1 = [eye(n); zeros(n,n)]; % selection matrix of first solution state T1 = [eye(n); zeros(n,n)]; % selection matrix of first solution state
T2 = [zeros(n,n); eye(n,n)]; % selection matrix of second solution state T2 = [zeros(n,n); eye(n,n)]; % selection matrix of second solution state
...@@ -193,7 +211,7 @@ for jBlock = 1:length(sigMod.blockLength) ...@@ -193,7 +211,7 @@ for jBlock = 1:length(sigMod.blockLength)
for zeta_idx=1:z_idx for zeta_idx=1:z_idx
temp_prod(zeta_idx,:,:) = ... temp_prod(zeta_idx,:,:) = ...
reshape(Psi_mu(z_idx,zeta_idx,:,:),[2*n 2*n])... reshape(Psi_mu(z_idx,zeta_idx,:,:),[2*n 2*n])...
*T2*(reshape(LDisc(zeta_idx,:,:),[n n])\reshape(A0(zeta_idx,:,:),[n n]))*T1.'; *T2*(reshape(LDisc(zeta_idx,:,:),[n n])\reshape(A0Disc(zeta_idx,:,:),[n n]))*T1.';
end end
A0_til(z_idx,:,:) = trapz(zdisc(1:z_idx),temp_prod,1); A0_til(z_idx,:,:) = trapz(zdisc(1:z_idx),temp_prod,1);
Y(z_idx,:,:) = reshape(Psi_mu(z_idx,1,:,:),[2*n,2*n]) - reshape(A0_til(z_idx,:,:),[2*n 2*n]); Y(z_idx,:,:) = reshape(Psi_mu(z_idx,1,:,:),[2*n,2*n]) - reshape(A0_til(z_idx,:,:),[2*n 2*n]);
...@@ -226,7 +244,7 @@ for jBlock = 1:length(sigMod.blockLength) ...@@ -226,7 +244,7 @@ for jBlock = 1:length(sigMod.blockLength)
g2jk = G2*varphi_j_k; g2jk = G2*varphi_j_k;
for z_idx = 1:ndisc for z_idx = 1:ndisc
hjk(z_idx,:) = reshape(H(z_idx,:,:),[n nv])*varphi_j_k; hjk(z_idx,:) = reshape(HDisc(z_idx,:,:),[n nv])*varphi_j_k;
temp_prod = zeros(z_idx,2*n,1); % preallocate and reset temp_prod = zeros(z_idx,2*n,1); % preallocate and reset
for zeta_idx=1:z_idx for zeta_idx=1:z_idx
temp_prod(zeta_idx,:,1) = ... temp_prod(zeta_idx,:,1) = ...
...@@ -323,10 +341,10 @@ for i=1:n ...@@ -323,10 +341,10 @@ for i=1:n
end end
end end
for z_idx=1:ndisc for z_idx=1:ndisc
lhs(z_idx,:,:) = sd(L(z_idx,:,:))*sd(Xzz(z_idx,:,:))... lhs(z_idx,:,:) = sd(LDisc(z_idx,:,:))*sd(Xzz(z_idx,:,:))...
+ sd(M(z_idx,:,:))*sd(sol(z_idx,:,:))... + sd(muSys(z_idx,:,:))*sd(sol(z_idx,:,:))...
+ sd(sol(z_idx,:,:))*S ... + sd(sol(z_idx,:,:))*S ...
+ sd(A0(z_idx,:,:))*sd(sol(1,:,:)); + sd(A0Disc(z_idx,:,:))*sd(sol(1,:,:));
end end
% rhs = H; % rhs = H;
% err = lhs-rhs; % err = lhs-rhs;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment