Mercurial > repos > public > sbplib

annotate +sbp/blocknorm4.m @ 132:6ec2248b83c4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactored +anim.animate to be more clear.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 05 Feb 2016 10:49:16 +0100
parents 32b39dc44474
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
29
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 function [H, HI, D1, D2, e_1, e_m, M,Q S_1, S_m] = blocknorm4(m,h)
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 H_U=[0.751e3 / 0.3456e4 0.661e3 / 0.3456e4 -0.515e3 / 0.3456e4 0.5e1 / 0.128e3; 0.661e3 / 0.3456e4 0.1405e4 / 0.1152e4 -0.3e1 / 0.128e3 0.29e2 / 0.3456e4; -0.515e3 / 0.3456e4 -0.3e1 / 0.128e3 0.989e3 / 0.1152e4 0.149e3 / 0.3456e4; 0.5e1 / 0.128e3 0.29e2 / 0.3456e4 0.149e3 / 0.3456e4 0.3407e4 / 0.3456e4;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 H=eye(m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 H(1:4,1:4)=H_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 H(m-3:m,m-3:m)=flipud( fliplr(H_U(1:4,1:4) ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 H=H*h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 HI=inv(H);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 Q_U = [-0.1e1 / 0.2e1 0.55e2 / 0.72e2 -0.47e2 / 0.144e3 0.1e1 / 0.16e2; -0.55e2 / 0.72e2 0 0.43e2 / 0.48e2 -0.19e2 / 0.144e3; 0.47e2 / 0.144e3 -0.43e2 / 0.48e2 0 0.47e2 / 0.72e2; -0.1e1 / 0.16e2 0.19e2 / 0.144e3 -0.47e2 / 0.72e2 0;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 Q(1:4,1:4)=Q_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 D1=HI*Q;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 M_U=[0.359e3 / 0.288e3 -0.443e3 / 0.288e3 0.97e2 / 0.288e3 -0.13e2 / 0.288e3; -0.51e2 / 0.32e2 0.325e3 / 0.96e2 -0.191e3 / 0.96e2 0.19e2 / 0.96e2; 0.43e2 / 0.96e2 -0.69e2 / 0.32e2 0.293e3 / 0.96e2 -0.137e3 / 0.96e2; -0.29e2 / 0.288e3 0.89e2 / 0.288e3 -0.427e3 / 0.288e3 0.727e3 / 0.288e3;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 M(1:4,1:4)=M_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 M(m-3:m,m-3:m)=flipud( fliplr( M_U ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29 M=M/h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 DS_U=[0.25e2 / 0.12e2 -4 3 -0.4e1 / 0.3e1 0.1e1 / 0.4e1;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 DS=zeros(m,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 DS(1,1:5)=DS_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 DS(m,m-4:m)=fliplr(DS_U);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 DS=DS/h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 D2=HI*(-M+DS);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 d3=[-1 3 -3 1];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 t3=sum(abs(d3));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 DD_3(1:1,1:4)=[d3];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 DD_3(m:m,m-3:m)=[d3];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 % This works for wave eq.
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 % For studs interface in 1D no AD is needed.
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 ADD=1*h/(t3)*DD_3'*DD_3;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 e_1 = zeros(m,1);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 e_1(1)= 1;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 e_m = zeros(m,1);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 e_m(end)= 1;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 S_1 = -DS(1,:)';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 S_m = DS(end,:)';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 Q = H*D1-(-e_1*e_1' + e_m*e_m');
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 end