Mercurial > repos > public > sbplib

annotate +sbp/higher_compatible6.m @ 255:df3cc9c5dffc operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added ordinary 12th order accurate, with D1*D1 as 2nd derivative.
author Martin Almquist <martin.almquist@it.uu.se>
date Wed, 07 Sep 2016 15:54:41 +0200
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, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = higher_compatible6(m,h)
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 %%% 4:de ordn. SBP Finita differens %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 %%% operatorer framtagna av Ken Mattsson %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 %%% %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 %%% 6 randpunkter, diagonal norm %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 %%% %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 %%% Datum: 2013-11-11 %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 %%% %%%
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 %%% H (Normen) %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 %%% D1 (approx f?rsta derivatan) %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 %%% D2 (approx andra derivatan) %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 %%% D3 (approx tredje derivatan) %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 %%% D2 (approx fj?rde derivatan) %%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 %%% %%%
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
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 % M?ste ange antal punkter (m) och stegl?ngd (h)
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 % Notera att dessa opetratorer ?r framtagna f?r att anv?ndas n?r
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 % vi har 3de och 4de derivator i v?r PDE
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 % I annat fall anv?nd de "traditionella" som har noggrannare
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 % randsplutningar f?r D1 och D2
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm,
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 % vilket ?r n?dv?ndigt f?r stabilitet
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 H=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
29 H_U=[0.7493827e7 / 0.25401600e8 0 0 0 0 0 0 0; 0 0.5534051e7 / 0.3628800e7 0 0 0 0 0 0; 0 0 0.104561e6 / 0.403200e6 0 0 0 0 0; 0 0 0 0.260503e6 / 0.145152e6 0 0 0 0; 0 0 0 0 0.43237e5 / 0.103680e6 0 0 0; 0 0 0 0 0 0.514081e6 / 0.403200e6 0 0; 0 0 0 0 0 0 0.3356179e7 / 0.3628800e7 0; 0 0 0 0 0 0 0 0.25631027e8 / 0.25401600e8;];
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 H(1:8,1:8)=H_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 H(m-7:m,m-7:m)=fliplr(flipud(H_U));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 H=H*h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 HI=inv(H);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35
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 % First derivative SBP operator, 3rd order accurate at first 8 boundary points
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 q3=1/60;q2=-3/20;q1=3/4;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2) - diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 Q_U = [0 0.26431903e8 / 0.43545600e8 0.39791489e8 / 0.152409600e9 -0.12747751e8 / 0.16934400e8 0.76099447e8 / 0.152409600e9 -0.1397443e7 / 0.12192768e8 0 0; -0.26431903e8 / 0.43545600e8 0 -0.13847476213e11 / 0.19559232000e11 0.35844843977e11 / 0.11735539200e11 -0.63413503537e11 / 0.23471078400e11 0.4764412871e10 / 0.3911846400e10 -0.1668252557e10 / 0.5867769600e10 0.842644697e9 / 0.29338848000e11; -0.39791489e8 / 0.152409600e9 0.13847476213e11 / 0.19559232000e11 0 -0.73834802771e11 / 0.23471078400e11 0.1802732209e10 / 0.325987200e9 -0.65514173e8 / 0.16299360e8 0.79341409141e11 / 0.58677696000e11 -0.1282384321e10 / 0.7823692800e10; 0.12747751e8 / 0.16934400e8 -0.35844843977e11 / 0.11735539200e11 0.73834802771e11 / 0.23471078400e11 0 -0.5274106087e10 / 0.1173553920e10 0.33743985841e11 / 0.5867769600e10 -0.6482602549e10 / 0.2607897600e10 0.1506017269e10 / 0.3911846400e10; -0.76099447e8 / 0.152409600e9 0.63413503537e11 / 0.23471078400e11 -0.1802732209e10 / 0.325987200e9 0.5274106087e10 / 0.1173553920e10 0 -0.7165829063e10 / 0.2607897600e10 0.23903110999e11 / 0.11735539200e11 -0.5346675911e10 / 0.11735539200e11; 0.1397443e7 / 0.12192768e8 -0.4764412871e10 / 0.3911846400e10 0.65514173e8 / 0.16299360e8 -0.33743985841e11 / 0.5867769600e10 0.7165829063e10 / 0.2607897600e10 0 -0.1060918223e10 / 0.11735539200e11 0.628353989e9 / 0.3911846400e10; 0 0.1668252557e10 / 0.5867769600e10 -0.79341409141e11 / 0.58677696000e11 0.6482602549e10 / 0.2607897600e10 -0.23903110999e11 / 0.11735539200e11 0.1060918223e10 / 0.11735539200e11 0 0.25889988599e11 / 0.39118464000e11; 0 -0.842644697e9 / 0.29338848000e11 0.1282384321e10 / 0.7823692800e10 -0.1506017269e10 / 0.3911846400e10 0.5346675911e10 / 0.11735539200e11 -0.628353989e9 / 0.3911846400e10 -0.25889988599e11 / 0.39118464000e11 0;];
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 Q(1:8,1:8)=Q_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 Q(m-7:m,m-7:m)=flipud( fliplr( -Q_U ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 e_1=zeros(m,1);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
48 e_m=zeros(m,1);e_m(m)=1;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 D1=HI*(Q-1/2*e_1*e_1'+1/2*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
52
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 % Second derivative, 1st order accurate at first 6 boundary points
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 % m3=-1/90;m2=3/20;m1=-3/2;m0=49/18;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 % M=m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*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
61 % M_U=[0.4347276223e10 / 0.3736212480e10 -0.1534657609e10 / 0.1210809600e10 0.68879e5 / 0.3057600e7 0.1092927401e10 / 0.13076743680e11 0.18145423e8 / 0.968647680e9 -0.1143817e7 / 0.60540480e8 -0.355447739e9 / 0.65383718400e11 0.56081e5 / 0.16473600e8; -0.1534657609e10 / 0.1210809600e10 0.42416226217e11 / 0.18681062400e11 -0.228654119e9 / 0.345945600e9 -0.12245627e8 / 0.34594560e8 -0.2995295e7 / 0.46702656e8 0.52836503e8 / 0.691891200e9 0.119351e6 / 0.12812800e8 -0.634102039e9 / 0.65383718400e11; 0.68879e5 / 0.3057600e7 -0.228654119e9 / 0.345945600e9 0.5399287e7 / 0.4193280e7 -0.24739409e8 / 0.34594560e8 0.7878737e7 / 0.69189120e8 -0.1917829e7 / 0.31449600e8 0.39727e5 / 0.3660800e7 0.10259e5 / 0.4656960e7; 0.1092927401e10 / 0.13076743680e11 -0.12245627e8 / 0.34594560e8 -0.24739409e8 / 0.34594560e8 0.7780367599e10 / 0.3736212480e10 -0.70085363e8 / 0.69189120e8 -0.500209e6 / 0.6289920e7 -0.311543e6 / 0.17962560e8 0.278191e6 / 0.21525504e8; 0.18145423e8 / 0.968647680e9 -0.2995295e7 / 0.46702656e8 0.7878737e7 / 0.69189120e8 -0.70085363e8 / 0.69189120e8 0.7116321131e10 / 0.3736212480e10 -0.545081e6 / 0.532224e6 0.811631e6 / 0.11531520e8 -0.84101639e8 / 0.13076743680e11; -0.1143817e7 / 0.60540480e8 0.52836503e8 / 0.691891200e9 -0.1917829e7 / 0.31449600e8 -0.500209e6 / 0.6289920e7 -0.545081e6 / 0.532224e6 0.324760747e9 / 0.138378240e9 -0.65995697e8 / 0.49420800e8 0.1469203e7 / 0.13759200e8; -0.355447739e9 / 0.65383718400e11 0.119351e6 / 0.12812800e8 0.39727e5 / 0.3660800e7 -0.311543e6 / 0.17962560e8 0.811631e6 / 0.11531520e8 -0.65995697e8 / 0.49420800e8 0.48284442317e11 / 0.18681062400e11 -0.1762877569e10 / 0.1210809600e10; 0.56081e5 / 0.16473600e8 -0.634102039e9 / 0.65383718400e11 0.10259e5 / 0.4656960e7 0.278191e6 / 0.21525504e8 -0.84101639e8 / 0.13076743680e11 0.1469203e7 / 0.13759200e8 -0.1762877569e10 / 0.1210809600e10 0.10117212851e11 / 0.3736212480e10;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 % M(1:8,1:8)=M_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 % M(m-7:m,m-7: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
66 % M=M/h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 S_U=[-0.12700800e8 / 0.7493827e7 0.185023321e9 / 0.89925924e8 0.39791489e8 / 0.44962962e8 -0.38243253e8 / 0.14987654e8 0.76099447e8 / 0.44962962e8 -0.34936075e8 / 0.89925924e8;]/h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 S_1=zeros(1,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 S_1(1:6)=S_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 S_m=zeros(1,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 S_m(m-5:m)=fliplr(-S_U);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75 %D2=HI*(-M-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
76
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 % Third derivative, 1st order accurate at first 6 boundary points
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84 % q4=7/240;q3=-3/10;q2=169/120;q1=-61/30;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 % Q3=q4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+q3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 % %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 % Q3_U = [0 -0.10882810591e11 / 0.5811886080e10 0.398713069e9 / 0.132088320e9 -0.1746657571e10 / 0.1162377216e10 0.56050639e8 / 0.145297152e9 -0.11473393e8 / 0.1162377216e10 -0.38062741e8 / 0.1452971520e10 0.30473e5 / 0.4392960e7; 0.10882810591e11 / 0.5811886080e10 0 -0.3720544343e10 / 0.830269440e9 0.767707019e9 / 0.207567360e9 -0.1047978301e10 / 0.830269440e9 0.1240729e7 / 0.14826240e8 0.6807397e7 / 0.55351296e8 -0.50022767e8 / 0.1452971520e10; -0.398713069e9 / 0.132088320e9 0.3720544343e10 / 0.830269440e9 0 -0.2870078009e10 / 0.830269440e9 0.74962049e8 / 0.29652480e8 -0.12944857e8 / 0.30750720e8 -0.17846623e8 / 0.103783680e9 0.68707591e8 / 0.1162377216e10; 0.1746657571e10 / 0.1162377216e10 -0.767707019e9 / 0.207567360e9 0.2870078009e10 / 0.830269440e9 0 -0.727867087e9 / 0.276756480e9 0.327603877e9 / 0.207567360e9 -0.175223717e9 / 0.830269440e9 0.1353613e7 / 0.726485760e9; -0.56050639e8 / 0.145297152e9 0.1047978301e10 / 0.830269440e9 -0.74962049e8 / 0.29652480e8 0.727867087e9 / 0.276756480e9 0 -0.1804641793e10 / 0.830269440e9 0.311038417e9 / 0.207567360e9 -0.1932566239e10 / 0.5811886080e10; 0.11473393e8 / 0.1162377216e10 -0.1240729e7 / 0.14826240e8 0.12944857e8 / 0.30750720e8 -0.327603877e9 / 0.207567360e9 0.1804641793e10 / 0.830269440e9 0 -0.1760949511e10 / 0.830269440e9 0.2105883973e10 / 0.1452971520e10; 0.38062741e8 / 0.1452971520e10 -0.6807397e7 / 0.55351296e8 0.17846623e8 / 0.103783680e9 0.175223717e9 / 0.830269440e9 -0.311038417e9 / 0.207567360e9 0.1760949511e10 / 0.830269440e9 0 -0.1081094773e10 / 0.528353280e9; -0.30473e5 / 0.4392960e7 0.50022767e8 / 0.1452971520e10 -0.68707591e8 / 0.1162377216e10 -0.1353613e7 / 0.726485760e9 0.1932566239e10 / 0.5811886080e10 -0.2105883973e10 / 0.1452971520e10 0.1081094773e10 / 0.528353280e9 0;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 % Q3(1:8,1:8)=Q3_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 % Q3(m-7:m,m-7:m)=flipud( fliplr( -Q3_U ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 % Q3=Q3/h^2;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 S2_U=[0.35e2 / 0.12e2 -0.26e2 / 0.3e1 0.19e2 / 0.2e1 -0.14e2 / 0.3e1 0.11e2 / 0.12e2;]/h^2;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 S2_1=zeros(1,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 S2_1(1:5)=S2_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 S2_m=zeros(1,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 S2_m(m-4:m)=fliplr(S2_U);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 %
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106 % D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 % Fourth derivative, 0th order accurate at first 6 boundary points (still
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 m4=7/240;m3=-2/5;m2=169/60;m1=-122/15;m0=91/8;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112 M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*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
113
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
115
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
116 M4_U=[0.600485868980522851e18 / 0.274825314114120000e18 -0.1421010223681841e16 / 0.348984525859200e15 0.38908412970187e14 / 0.1586293299360000e16 0.10224077451922837e17 / 0.2243471951952000e16 -0.7577302712815639e16 / 0.1744922629296000e16 0.138091642084013e15 / 0.59351109840000e14 -0.3775041725375197e16 / 0.4486943903904000e16 0.9907210230881393e16 / 0.61072292025360000e17; -0.1421010223681841e16 / 0.348984525859200e15 0.3985852497808703e16 / 0.407903991264000e15 -0.90048788923861e14 / 0.15579666333000e14 -0.4312795866499e13 / 0.997098645312e12 0.4414634708891947e16 / 0.448694390390400e15 -0.886174803100459e15 / 0.99709864531200e14 0.4333e4 / 0.1000e4 -0.13800578064893047e17 / 0.15704303663664000e17; 0.38908412970187e14 / 0.1586293299360000e16 -0.90048788923861e14 / 0.15579666333000e14 0.2071682582321887e16 / 0.113306664240000e15 -0.769471337294003e15 / 0.41545776888000e14 0.112191585452033e15 / 0.166183107552000e15 0.7204491902193671e16 / 0.623186653320000e15 -0.24847093554379e14 / 0.3115933266600e13 0.943854037768721e15 / 0.545288321655000e15; 0.10224077451922837e17 / 0.2243471951952000e16 -0.4312795866499e13 / 0.997098645312e12 -0.769471337294003e15 / 0.41545776888000e14 0.3086874339649421e16 / 0.81580798252800e14 -0.396009005312111e15 / 0.16618310755200e14 0.348854811893087e15 / 0.249274661328000e15 0.895954627955053e15 / 0.224347195195200e15 -0.184881685054543e15 / 0.166183107552000e15; -0.7577302712815639e16 / 0.1744922629296000e16 0.4414634708891947e16 / 0.448694390390400e15 0.112191585452033e15 / 0.166183107552000e15 -0.396009005312111e15 / 0.16618310755200e14 0.3774861828677557e16 / 0.112173597597600e15 -0.5693689108983593e16 / 0.249274661328000e15 0.803944126167107e15 / 0.99709864531200e14 -0.19547569411550791e17 / 0.15704303663664000e17; 0.138091642084013e15 / 0.59351109840000e14 -0.886174803100459e15 / 0.99709864531200e14 0.7204491902193671e16 / 0.623186653320000e15 0.348854811893087e15 / 0.249274661328000e15 -0.5693689108983593e16 / 0.249274661328000e15 0.73965842628398389e17 / 0.2492746613280000e16 -0.2184472662036043e16 / 0.124637330664000e15 0.46667e5 / 0.10000e5; -0.3775041725375197e16 / 0.4486943903904000e16 0.4333e4 / 0.1000e4 -0.24847093554379e14 / 0.3115933266600e13 0.895954627955053e15 / 0.224347195195200e15 0.803944126167107e15 / 0.99709864531200e14 -0.2184472662036043e16 / 0.124637330664000e15 0.37593640125444199e17 / 0.2243471951952000e16 -0.37e2 / 0.4e1; 0.9907210230881393e16 / 0.61072292025360000e17 -0.13800578064893047e17 / 0.15704303663664000e17 0.943854037768721e15 / 0.545288321655000e15 -0.184881685054543e15 / 0.166183107552000e15 -0.19547569411550791e17 / 0.15704303663664000e17 0.46667e5 / 0.10000e5 -0.37e2 / 0.4e1 0.12766926490502478779e20 / 0.1099301256456480000e19;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
117
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
118 M4(1:8,1:8)=M4_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
119
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
120 M4(m-7:m,m-7:m)=flipud( fliplr( M4_U ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
121 M4=M4/h^3;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
122
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
123 S3_U=[-0.5e1 / 0.2e1 9 -12 7 -0.3e1 / 0.2e1;]/h^3;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
124 S3_1=zeros(1,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
125 S3_1(1:5)=S3_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
126 S3_m=zeros(1,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
127 S3_m(m-4:m)=fliplr(-S3_U);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
128
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
129 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
130
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
131 S_1 = S_1';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
132 S_m = S_m';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
133 S2_1 = S2_1';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
134 S2_m = S2_m';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
135 S3_1 = S3_1';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
136 S3_m = S3_m';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
137
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
138 % L=h*(m-1);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
139
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
140 % x1=linspace(0,L,m)';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
141 % x2=x1.^2/fac(2);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
142 % x3=x1.^3/fac(3);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
143 % x4=x1.^4/fac(4);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
144 % x5=x1.^5/fac(5);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
145
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
146 % x0=x1.^0/fac(1);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
147
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
148 end