Mercurial > repos > public > sbplib

comparison +sbp/ordinary2.m @ 29:32b39dc44474

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Removed repository inside +sbp to make it part of the root repo.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 28 Sep 2015 08:47:28 +0200
parents
children
comparison
equal deleted inserted replaced
28:16acb2775aca 29:32b39dc44474
1 function [H, HI, D1, D2, e_1, e_m, M,Q S_1, S_m] = ordinary2(m,h)
2 e_1=zeros(m,1);e_1(1)=1;
3 e_m=zeros(m,1);e_m(m)=1;
4
5 H=(eye(m,m));H(1,1)=0.5;H(m,m)=0.5;
6 H=h*H;
7 HI=inv(H);
8
9 D1=((.5*diag(ones(m-1,1),1)-.5*diag(ones(m-1,1),-1)));
10 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1;
11 D1(m,m-1)=-1;D1(m,m)=1;
12 D1=D1/h;
13
14 Q=H*D1 + 1/2*e_1*e_1' - 1/2*e_m*e_m';
15
16 D2=((diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)-2*diag(ones(m,1),0)));
17 D2(1,1)=1;D2(1,2)=-2;D2(1,3)=1;
18 D2(m,m-2)=1;D2(m,m-1)=-2;D2(m,m)=1;
19 D2=D2/h^2;
20
21 S_U=[-3/2, 2, -1/2]/h;
22 S_1=zeros(1,m);
23 S_1(1:3)=S_U;
24 S_m=zeros(1,m);
25 S_m(m-2:m)=fliplr(-S_U);
26
27
28 M=-H*D2-e_1*S_1+e_m*S_m;
29 S_1 = S_1';
30 S_m = S_m';
31 end