Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_4.m @ 315:297d2cbfbe15 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleaning of variable 2nd deriv in 4th order.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 22:58:45 +0200
parents 88584b0cfba1
children b19e142fcae1
comparison
equal deleted inserted replaced
314:88584b0cfba1 315:297d2cbfbe15
57 57
58 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); 58 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r');
59 59
60 60
61 % Second derivative 61 % Second derivative
62 M=sparse(m,m); 62 M = sparse(m,m);
63 63
64 scheme_width = 5; 64 scheme_width = 5;
65 scheme_radius = (scheme_width-1)/2; 65 scheme_radius = (scheme_width-1)/2;
66 r = (1+scheme_radius):(m-scheme_radius); 66 r = (1+scheme_radius):(m-scheme_radius);
67 67
68 function D2 = D2_fun(c) 68 function D2 = D2_fun(c)
69 % ALTERNATIVES %%%%%%%%%%%%%
70 % for i=4:m-3
71 % M(i,i-2:i+2)=[-c(i-1)/0.6e1 + c(i-2)/0.8e1 + c(i)/0.8e1 -c(i-2)/0.6e1 - c(i+1)/0.6e1 - c(i-1)/0.2e1 - c(i)/0.2e1 c(i-2)/0.24e2 + 0.5e1/0.6e1 * c(i-1) + 0.5e1/0.6e1 * c(i+1) + c(i+2)/0.24e2 + 0.3e1/0.4e1 * c(i) -c(i-1)/0.6e1 - c(i+2)/0.6e1 - c(i)/0.2e1 - c(i+1)/0.2e1 -c(i+1)/0.6e1 + c(i)/0.8e1 + c(i+2)/0.8e1;];
72 % end
73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
74 % for i=4:m-3
75 % M(i,i-2:i+2)= [
76 % c(i-2)/0.8e1 - c(i-1)/0.6e1 + c(i) /0.8e1 ,
77 % -c(i-2)/0.6e1 - c(i-1)/0.2e1 - c(i) /0.2e1 - c(i+1)/0.6e1 ,
78 % c(i-2)/2.4e1 + c(i-1)/1.2e0 + c(i) * 0.3/0.4 + c(i+1)/1.2e0 + c(i+2)/2.4e1 ,
79 % -c(i-1)/0.6e1 - c(i) /0.2e1 - c(i+1)/0.2e1 - c(i+2)/0.6e1 ,
80 % c(i) /0.8e1 - c(i+1)/0.6e1 + c(i+2)/0.8e1 ,
81 % ];
82 % end
83 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84 Mm2 = c(r-2)/0.8e1 - c(r-1)/0.6e1 + c(r) /0.8e1 ; 69 Mm2 = c(r-2)/0.8e1 - c(r-1)/0.6e1 + c(r) /0.8e1 ;
85 Mm1 = -c(r-2)/0.6e1 - c(r-1)/0.2e1 - c(r) /0.2e1 - c(r+1)/0.6e1 ; 70 Mm1 = -c(r-2)/0.6e1 - c(r-1)/0.2e1 - c(r) /0.2e1 - c(r+1)/0.6e1 ;
86 M0 = c(r-2)/2.4e1 + c(r-1)/1.2e0 + c(r) * 0.3/0.4 + c(r+1)/1.2e0 + c(r+2)/2.4e1; 71 M0 = c(r-2)/2.4e1 + c(r-1)/1.2e0 + c(r) * 0.3/0.4 + c(r+1)/1.2e0 + c(r+2)/2.4e1;
87 Mp1 = -c(r-1)/0.6e1 - c(r) /0.2e1 - c(r+1)/0.2e1 - c(r+2)/0.6e1; 72 Mp1 = -c(r-1)/0.6e1 - c(r) /0.2e1 - c(r+1)/0.2e1 - c(r+2)/0.6e1;
88 Mp2 = c(r) /0.8e1 - c(r+1)/0.6e1 + c(r+2)/0.8e1; 73 Mp2 = c(r) /0.8e1 - c(r+1)/0.6e1 + c(r+2)/0.8e1;
89 % printSize(Mm2); 74
90 % scheme_radius
91 % m
92 M(r,:) = spdiags([Mm2 Mm1 M0 Mp1 Mp2],0:2*scheme_radius,length(r),m); 75 M(r,:) = spdiags([Mm2 Mm1 M0 Mp1 Mp2],0:2*scheme_radius,length(r),m);
93 % M(r,:) = spdiags([Mm2 Mm1 M0 Mp1 Mp2],(-2:2)+scheme_radius,M(r,:)); % This is slower
94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95 % %% Somthing is wrong here!!
96 % Mm2 = c(r-2)/0.8e1 - c(r-1)/0.6e1 + c(r) /0.8e1 ;
97 % Mm1 = -c(r-2)/0.6e1 - c(r-1)/0.2e1 - c(r) /0.2e1 - c(r+1)/0.6e1 ;
98 % M0 = c(r-2)/2.4e1 + c(r-1)/1.2e0 + c(r) * 0.3/0.4 + c(r+1)/1.2e0 + c(r+2)/2.4e1;
99 % Mp1 = -c(r-1)/0.6e1 - c(r) /0.2e1 - c(r+1)/0.2e1 - c(r+2)/0.6e1;
100 % Mp2 = c(r) /0.8e1 - c(r+1)/0.6e1 + c(r+2)/0.8e1;
101 % % printSize(M_diag_ind);
102 % % Mdiags = [Mm2 Mm1 M0 Mp1 Mp2];
103 % % printSize(Mdiags);
104 % M(M_diag_ind) = [Mm2 Mm1 M0 Mp1 Mp2]; % This is slightly faster
105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 % Kan man skriva det som en multiplikation av en 3-dim matris?
107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108 76
109 M(1:6,1:6) = [ 77 M(1:6,1:6) = [
110 0.12e2/0.17e2 * c(1) + 0.59e2/0.192e3 * c(2) + 0.27010400129e11/0.345067064608e12 * c(3) + 0.69462376031e11/0.2070402387648e13 * c(4) -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -c(4)/0.784e3 + c(3)/0.784e3; 78 0.12e2/0.17e2 * c(1) + 0.59e2/0.192e3 * c(2) + 0.27010400129e11/0.345067064608e12 * c(3) + 0.69462376031e11/0.2070402387648e13 * c(4) -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -c(4)/0.784e3 + c(3)/0.784e3;
111 -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.3481e4/0.3264e4 * c(1) + 0.9258282831623875e16/0.7669235228057664e16 * c(3) + 0.236024329996203e15/0.1278205871342944e16 * c(4) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4); 79 -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.3481e4/0.3264e4 * c(1) + 0.9258282831623875e16/0.7669235228057664e16 * c(3) + 0.236024329996203e15/0.1278205871342944e16 * c(4) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4);
112 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) c(1)/0.51e2 + 0.59e2/0.192e3 * c(2) + 0.13777050223300597e17/0.26218083221499456e17 * c(4) + 0.564461e6/0.13384296e8 * c(5) + 0.378288882302546512209e21/0.270764341349677687456e21 * c(3) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) -0.10532412077335e14/0.42840005263888e14 * c(4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(3) + 0.564461e6/0.4461432e7 * c(5) -0.960119e6/0.1280713392e10 * c(4) - 0.3391e4/0.6692148e7 * c(5) + 0.33235054191e11/0.26452850508784e14 * c(3); 80 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) c(1)/0.51e2 + 0.59e2/0.192e3 * c(2) + 0.13777050223300597e17/0.26218083221499456e17 * c(4) + 0.564461e6/0.13384296e8 * c(5) + 0.378288882302546512209e21/0.270764341349677687456e21 * c(3) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) -0.10532412077335e14/0.42840005263888e14 * c(4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(3) + 0.564461e6/0.4461432e7 * c(5) -0.960119e6/0.1280713392e10 * c(4) - 0.3391e4/0.6692148e7 * c(5) + 0.33235054191e11/0.26452850508784e14 * c(3);
113 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) 0.3e1/0.1088e4 * c(1) + 0.507284006600757858213e21/0.475219048083107777984e21 * c(3) + 0.1869103e7/0.2230716e7 * c(5) + c(6)/0.24e2 + 0.1950062198436997e16/0.3834617614028832e16 * c(4) -0.4959271814984644613e19/0.20965546238960637264e20 * c(3) - c(6)/0.6e1 - 0.15998714909649e14/0.37594290333616e14 * c(4) - 0.375177e6/0.743572e6 * c(5) -0.368395e6/0.2230716e7 * c(5) + 0.752806667e9/0.539854092016e12 * c(3) + 0.1063649e7/0.8712336e7 * c(4) + c(6)/0.8e1; 81 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) 0.3e1/0.1088e4 * c(1) + 0.507284006600757858213e21/0.475219048083107777984e21 * c(3) + 0.1869103e7/0.2230716e7 * c(5) + c(6)/0.24e2 + 0.1950062198436997e16/0.3834617614028832e16 * c(4) -0.4959271814984644613e19/0.20965546238960637264e20 * c(3) - c(6)/0.6e1 - 0.15998714909649e14/0.37594290333616e14 * c(4) - 0.375177e6/0.743572e6 * c(5) -0.368395e6/0.2230716e7 * c(5) + 0.752806667e9/0.539854092016e12 * c(3) + 0.1063649e7/0.8712336e7 * c(4) + c(6)/0.8e1;
114 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.10532412077335e14/0.42840005263888e14 * c(4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(3) + 0.564461e6/0.4461432e7 * c(5) -0.4959271814984644613e19/0.20965546238960637264e20 * c(3) - c(6)/0.6e1 - 0.15998714909649e14/0.37594290333616e14 * c(4) - 0.375177e6/0.743572e6 * c(5) 0.8386761355510099813e19/0.128413970713633903242e21 * c(3) + 0.2224717261773437e16/0.2763180339520776e16 * c(4) + 0.5e1/0.6e1 * c(6) + c(7)/0.24e2 + 0.280535e6/0.371786e6 * c(5) -0.35039615e8/0.213452232e9 * c(4) - c(7)/0.6e1 - 0.13091810925e11/0.13226425254392e14 * c(3) - 0.1118749e7/0.2230716e7 * c(5) - c(6)/0.2e1; 82 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.10532412077335e14/0.42840005263888e14 * c(4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(3) + 0.564461e6/0.4461432e7 * c(5) -0.4959271814984644613e19/0.20965546238960637264e20 * c(3) - c(6)/0.6e1 - 0.15998714909649e14/0.37594290333616e14 * c(4) - 0.375177e6/0.743572e6 * c(5) 0.8386761355510099813e19/0.128413970713633903242e21 * c(3) + 0.2224717261773437e16/0.2763180339520776e16 * c(4) + 0.5e1/0.6e1 * c(6) + c(7)/0.24e2 + 0.280535e6/0.371786e6 * c(5) -0.35039615e8/0.213452232e9 * c(4) - c(7)/0.6e1 - 0.13091810925e11/0.13226425254392e14 * c(3) - 0.1118749e7/0.2230716e7 * c(5) - c(6)/0.2e1;
115 -c(4)/0.784e3 + c(3)/0.784e3 -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4) -0.960119e6/0.1280713392e10 * c(4) - 0.3391e4/0.6692148e7 * c(5) + 0.33235054191e11/0.26452850508784e14 * c(3) -0.368395e6/0.2230716e7 * c(5) + 0.752806667e9/0.539854092016e12 * c(3) + 0.1063649e7/0.8712336e7 * c(4) + c(6)/0.8e1 -0.35039615e8/0.213452232e9 * c(4) - c(7)/0.6e1 - 0.13091810925e11/0.13226425254392e14 * c(3) - 0.1118749e7/0.2230716e7 * c(5) - c(6)/0.2e1 0.3290636e7/0.80044587e8 * c(4) + 0.5580181e7/0.6692148e7 * c(5) + 0.5e1/0.6e1 * c(7) + c(8)/0.24e2 + 0.660204843e9/0.13226425254392e14 * c(3) + 0.3e1/0.4e1 * c(6); 83 -c(4)/0.784e3 + c(3)/0.784e3 -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4) -0.960119e6/0.1280713392e10 * c(4) - 0.3391e4/0.6692148e7 * c(5) + 0.33235054191e11/0.26452850508784e14 * c(3) -0.368395e6/0.2230716e7 * c(5) + 0.752806667e9/0.539854092016e12 * c(3) + 0.1063649e7/0.8712336e7 * c(4) + c(6)/0.8e1 -0.35039615e8/0.213452232e9 * c(4) - c(7)/0.6e1 - 0.13091810925e11/0.13226425254392e14 * c(3) - 0.1118749e7/0.2230716e7 * c(5) - c(6)/0.2e1 0.3290636e7/0.80044587e8 * c(4) + 0.5580181e7/0.6692148e7 * c(5) + 0.5e1/0.6e1 * c(7) + c(8)/0.24e2 + 0.660204843e9/0.13226425254392e14 * c(3) + 0.3e1/0.4e1 * c(6);
116 ]; 84 ];
117
118 85
119 M(m-5:m,m-5:m) = [ 86 M(m-5:m,m-5:m) = [
120 c(m-7)/0.24e2 + 0.5e1/0.6e1 * c(m-6) + 0.5580181e7/0.6692148e7 * c(m-4) + 0.4887707739997e13/0.119037827289528e15 * c(m-3) + 0.3e1/0.4e1 * c(m-5) + 0.660204843e9/0.13226425254392e14 * c(m-2) + 0.660204843e9/0.13226425254392e14 * c(m-1) -c(m-6)/0.6e1 - 0.1618585929605e13/0.9919818940794e13 * c(m-3) - c(m-5)/0.2e1 - 0.1118749e7/0.2230716e7 * c(m-4) - 0.13091810925e11/0.13226425254392e14 * c(m-2) - 0.13091810925e11/0.13226425254392e14 * c(m-1) -0.368395e6/0.2230716e7 * c(m-4) + c(m-5)/0.8e1 + 0.48866620889e11/0.404890569012e12 * c(m-3) + 0.752806667e9/0.539854092016e12 * c(m-2) + 0.752806667e9/0.539854092016e12 * c(m-1) -0.3391e4/0.6692148e7 * c(m-4) - 0.238797444493e12/0.119037827289528e15 * c(m-3) + 0.33235054191e11/0.26452850508784e14 * c(m-2) + 0.33235054191e11/0.26452850508784e14 * c(m-1) -0.8673e4/0.2904112e7 * c(m-2) - 0.8673e4/0.2904112e7 * c(m-1) + 0.8673e4/0.1452056e7 * c(m-3) -c(m-3)/0.392e3 + c(m-2)/0.784e3 + c(m-1)/0.784e3; 87 c(m-7)/0.24e2 + 0.5e1/0.6e1 * c(m-6) + 0.5580181e7/0.6692148e7 * c(m-4) + 0.4887707739997e13/0.119037827289528e15 * c(m-3) + 0.3e1/0.4e1 * c(m-5) + 0.660204843e9/0.13226425254392e14 * c(m-2) + 0.660204843e9/0.13226425254392e14 * c(m-1) -c(m-6)/0.6e1 - 0.1618585929605e13/0.9919818940794e13 * c(m-3) - c(m-5)/0.2e1 - 0.1118749e7/0.2230716e7 * c(m-4) - 0.13091810925e11/0.13226425254392e14 * c(m-2) - 0.13091810925e11/0.13226425254392e14 * c(m-1) -0.368395e6/0.2230716e7 * c(m-4) + c(m-5)/0.8e1 + 0.48866620889e11/0.404890569012e12 * c(m-3) + 0.752806667e9/0.539854092016e12 * c(m-2) + 0.752806667e9/0.539854092016e12 * c(m-1) -0.3391e4/0.6692148e7 * c(m-4) - 0.238797444493e12/0.119037827289528e15 * c(m-3) + 0.33235054191e11/0.26452850508784e14 * c(m-2) + 0.33235054191e11/0.26452850508784e14 * c(m-1) -0.8673e4/0.2904112e7 * c(m-2) - 0.8673e4/0.2904112e7 * c(m-1) + 0.8673e4/0.1452056e7 * c(m-3) -c(m-3)/0.392e3 + c(m-2)/0.784e3 + c(m-1)/0.784e3;
121 -c(m-6)/0.6e1 - 0.1618585929605e13/0.9919818940794e13 * c(m-3) - c(m-5)/0.2e1 - 0.1118749e7/0.2230716e7 * c(m-4) - 0.13091810925e11/0.13226425254392e14 * c(m-2) - 0.13091810925e11/0.13226425254392e14 * c(m-1) c(m-6)/0.24e2 + 0.5e1/0.6e1 * c(m-5) + 0.3896014498639e13/0.4959909470397e13 * c(m-3) + 0.8386761355510099813e19/0.128413970713633903242e21 * c(m-2) + 0.280535e6/0.371786e6 * c(m-4) + 0.3360696339136261875e19/0.171218627618178537656e21 * c(m-1) -c(m-5)/0.6e1 - 0.4959271814984644613e19/0.20965546238960637264e20 * c(m-2) - 0.375177e6/0.743572e6 * c(m-4) - 0.13425842714e11/0.33740880751e11 * c(m-3) - 0.193247108773400725e18/0.6988515412986879088e19 * c(m-1) -0.365281640980e12/0.1653303156799e13 * c(m-3) + 0.564461e6/0.4461432e7 * c(m-4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18/0.7963657098519931984e19 * c(m-1) -0.1328188692663e13/0.37594290333616e14 * c(m-2) + 0.2226377963775e13/0.37594290333616e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) c(m-3)/0.49e2 + 0.49579087e8/0.10149031312e11 * c(m-2) - 0.256702175e9/0.10149031312e11 * c(m-1); 88 -c(m-6)/0.6e1 - 0.1618585929605e13/0.9919818940794e13 * c(m-3) - c(m-5)/0.2e1 - 0.1118749e7/0.2230716e7 * c(m-4) - 0.13091810925e11/0.13226425254392e14 * c(m-2) - 0.13091810925e11/0.13226425254392e14 * c(m-1) c(m-6)/0.24e2 + 0.5e1/0.6e1 * c(m-5) + 0.3896014498639e13/0.4959909470397e13 * c(m-3) + 0.8386761355510099813e19/0.128413970713633903242e21 * c(m-2) + 0.280535e6/0.371786e6 * c(m-4) + 0.3360696339136261875e19/0.171218627618178537656e21 * c(m-1) -c(m-5)/0.6e1 - 0.4959271814984644613e19/0.20965546238960637264e20 * c(m-2) - 0.375177e6/0.743572e6 * c(m-4) - 0.13425842714e11/0.33740880751e11 * c(m-3) - 0.193247108773400725e18/0.6988515412986879088e19 * c(m-1) -0.365281640980e12/0.1653303156799e13 * c(m-3) + 0.564461e6/0.4461432e7 * c(m-4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18/0.7963657098519931984e19 * c(m-1) -0.1328188692663e13/0.37594290333616e14 * c(m-2) + 0.2226377963775e13/0.37594290333616e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) c(m-3)/0.49e2 + 0.49579087e8/0.10149031312e11 * c(m-2) - 0.256702175e9/0.10149031312e11 * c(m-1);
122 -0.368395e6/0.2230716e7 * c(m-4) + c(m-5)/0.8e1 + 0.48866620889e11/0.404890569012e12 * c(m-3) + 0.752806667e9/0.539854092016e12 * c(m-2) + 0.752806667e9/0.539854092016e12 * c(m-1) -c(m-5)/0.6e1 - 0.4959271814984644613e19/0.20965546238960637264e20 * c(m-2) - 0.375177e6/0.743572e6 * c(m-4) - 0.13425842714e11/0.33740880751e11 * c(m-3) - 0.193247108773400725e18/0.6988515412986879088e19 * c(m-1) c(m-5)/0.24e2 + 0.1869103e7/0.2230716e7 * c(m-4) + 0.507284006600757858213e21/0.475219048083107777984e21 * c(m-2) + 0.3e1/0.1088e4 * c(m) + 0.31688435395e11/0.67481761502e11 * c(m-3) + 0.27769176016102795561e20/0.712828572124661666976e21 * c(m-1) -0.125059e6/0.743572e6 * c(m-4) + c(m)/0.136e3 - 0.23099342648e11/0.101222642253e12 * c(m-3) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18/0.5525802884687299744e19 * c(m-1) 0.260297319232891e15/0.2556411742685888e16 * c(m-2) - 0.59e2/0.1088e4 * c(m) - 0.106641839640553e15/0.1278205871342944e16 * c(m-1) + 0.26019e5/0.726028e6 * c(m-3) -0.1244724001e10/0.21126554976e11 * c(m-2) + 0.3e1/0.68e2 * c(m) + 0.752806667e9/0.21126554976e11 * c(m-1); 89 -0.368395e6/0.2230716e7 * c(m-4) + c(m-5)/0.8e1 + 0.48866620889e11/0.404890569012e12 * c(m-3) + 0.752806667e9/0.539854092016e12 * c(m-2) + 0.752806667e9/0.539854092016e12 * c(m-1) -c(m-5)/0.6e1 - 0.4959271814984644613e19/0.20965546238960637264e20 * c(m-2) - 0.375177e6/0.743572e6 * c(m-4) - 0.13425842714e11/0.33740880751e11 * c(m-3) - 0.193247108773400725e18/0.6988515412986879088e19 * c(m-1) c(m-5)/0.24e2 + 0.1869103e7/0.2230716e7 * c(m-4) + 0.507284006600757858213e21/0.475219048083107777984e21 * c(m-2) + 0.3e1/0.1088e4 * c(m) + 0.31688435395e11/0.67481761502e11 * c(m-3) + 0.27769176016102795561e20/0.712828572124661666976e21 * c(m-1) -0.125059e6/0.743572e6 * c(m-4) + c(m)/0.136e3 - 0.23099342648e11/0.101222642253e12 * c(m-3) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18/0.5525802884687299744e19 * c(m-1) 0.260297319232891e15/0.2556411742685888e16 * c(m-2) - 0.59e2/0.1088e4 * c(m) - 0.106641839640553e15/0.1278205871342944e16 * c(m-1) + 0.26019e5/0.726028e6 * c(m-3) -0.1244724001e10/0.21126554976e11 * c(m-2) + 0.3e1/0.68e2 * c(m) + 0.752806667e9/0.21126554976e11 * c(m-1);
123 -0.3391e4/0.6692148e7 * c(m-4) - 0.238797444493e12/0.119037827289528e15 * c(m-3) + 0.33235054191e11/0.26452850508784e14 * c(m-2) + 0.33235054191e11/0.26452850508784e14 * c(m-1) -0.365281640980e12/0.1653303156799e13 * c(m-3) + 0.564461e6/0.4461432e7 * c(m-4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18/0.7963657098519931984e19 * c(m-1) -0.125059e6/0.743572e6 * c(m-4) + c(m)/0.136e3 - 0.23099342648e11/0.101222642253e12 * c(m-3) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18/0.5525802884687299744e19 * c(m-1) 0.564461e6/0.13384296e8 * c(m-4) + 0.470299699916357e15/0.952302618316224e15 * c(m-3) + 0.550597048646198778781e21/0.1624586048098066124736e22 * c(m-1) + c(m)/0.51e2 + 0.378288882302546512209e21/0.270764341349677687456e21 * c(m-2) -0.59e2/0.408e3 * c(m) - 0.29294615794607e14/0.29725717938208e14 * c(m-2) - 0.2234477713167e13/0.29725717938208e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) -0.59e2/0.3136e4 * c(m-3) - 0.13249937023e11/0.48148892736e11 * c(m-1) + 0.2e1/0.17e2 * c(m) + 0.2083938599e10/0.8024815456e10 * c(m-2); 90 -0.3391e4/0.6692148e7 * c(m-4) - 0.238797444493e12/0.119037827289528e15 * c(m-3) + 0.33235054191e11/0.26452850508784e14 * c(m-2) + 0.33235054191e11/0.26452850508784e14 * c(m-1) -0.365281640980e12/0.1653303156799e13 * c(m-3) + 0.564461e6/0.4461432e7 * c(m-4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18/0.7963657098519931984e19 * c(m-1) -0.125059e6/0.743572e6 * c(m-4) + c(m)/0.136e3 - 0.23099342648e11/0.101222642253e12 * c(m-3) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18/0.5525802884687299744e19 * c(m-1) 0.564461e6/0.13384296e8 * c(m-4) + 0.470299699916357e15/0.952302618316224e15 * c(m-3) + 0.550597048646198778781e21/0.1624586048098066124736e22 * c(m-1) + c(m)/0.51e2 + 0.378288882302546512209e21/0.270764341349677687456e21 * c(m-2) -0.59e2/0.408e3 * c(m) - 0.29294615794607e14/0.29725717938208e14 * c(m-2) - 0.2234477713167e13/0.29725717938208e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) -0.59e2/0.3136e4 * c(m-3) - 0.13249937023e11/0.48148892736e11 * c(m-1) + 0.2e1/0.17e2 * c(m) + 0.2083938599e10/0.8024815456e10 * c(m-2);
124 -0.8673e4/0.2904112e7 * c(m-2) - 0.8673e4/0.2904112e7 * c(m-1) + 0.8673e4/0.1452056e7 * c(m-3) -0.1328188692663e13/0.37594290333616e14 * c(m-2) + 0.2226377963775e13/0.37594290333616e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) 0.260297319232891e15/0.2556411742685888e16 * c(m-2) - 0.59e2/0.1088e4 * c(m) - 0.106641839640553e15/0.1278205871342944e16 * c(m-1) + 0.26019e5/0.726028e6 * c(m-3) -0.59e2/0.408e3 * c(m) - 0.29294615794607e14/0.29725717938208e14 * c(m-2) - 0.2234477713167e13/0.29725717938208e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) 0.9258282831623875e16/0.7669235228057664e16 * c(m-2) + 0.3481e4/0.3264e4 * c(m) + 0.228389721191751e15/0.1278205871342944e16 * c(m-1) + 0.8673e4/0.1452056e7 * c(m-3) -0.6025413881e10/0.21126554976e11 * c(m-2) - 0.59e2/0.68e2 * c(m) - 0.537416663e9/0.7042184992e10 * c(m-1); 91 -0.8673e4/0.2904112e7 * c(m-2) - 0.8673e4/0.2904112e7 * c(m-1) + 0.8673e4/0.1452056e7 * c(m-3) -0.1328188692663e13/0.37594290333616e14 * c(m-2) + 0.2226377963775e13/0.37594290333616e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) 0.260297319232891e15/0.2556411742685888e16 * c(m-2) - 0.59e2/0.1088e4 * c(m) - 0.106641839640553e15/0.1278205871342944e16 * c(m-1) + 0.26019e5/0.726028e6 * c(m-3) -0.59e2/0.408e3 * c(m) - 0.29294615794607e14/0.29725717938208e14 * c(m-2) - 0.2234477713167e13/0.29725717938208e14 * c(m-1) - 0.8673e4/0.363014e6 * c(m-3) 0.9258282831623875e16/0.7669235228057664e16 * c(m-2) + 0.3481e4/0.3264e4 * c(m) + 0.228389721191751e15/0.1278205871342944e16 * c(m-1) + 0.8673e4/0.1452056e7 * c(m-3) -0.6025413881e10/0.21126554976e11 * c(m-2) - 0.59e2/0.68e2 * c(m) - 0.537416663e9/0.7042184992e10 * c(m-1);
125 -c(m-3)/0.392e3 + c(m-2)/0.784e3 + c(m-1)/0.784e3 c(m-3)/0.49e2 + 0.49579087e8/0.10149031312e11 * c(m-2) - 0.256702175e9/0.10149031312e11 * c(m-1) -0.1244724001e10/0.21126554976e11 * c(m-2) + 0.3e1/0.68e2 * c(m) + 0.752806667e9/0.21126554976e11 * c(m-1) -0.59e2/0.3136e4 * c(m-3) - 0.13249937023e11/0.48148892736e11 * c(m-1) + 0.2e1/0.17e2 * c(m) + 0.2083938599e10/0.8024815456e10 * c(m-2) -0.6025413881e10/0.21126554976e11 * c(m-2) - 0.59e2/0.68e2 * c(m) - 0.537416663e9/0.7042184992e10 * c(m-1) 0.3e1/0.3136e4 * c(m-3) + 0.27010400129e11/0.345067064608e12 * c(m-2) + 0.234566387291e12/0.690134129216e12 * c(m-1) + 0.12e2/0.17e2 * c(m); 92 -c(m-3)/0.392e3 + c(m-2)/0.784e3 + c(m-1)/0.784e3 c(m-3)/0.49e2 + 0.49579087e8/0.10149031312e11 * c(m-2) - 0.256702175e9/0.10149031312e11 * c(m-1) -0.1244724001e10/0.21126554976e11 * c(m-2) + 0.3e1/0.68e2 * c(m) + 0.752806667e9/0.21126554976e11 * c(m-1) -0.59e2/0.3136e4 * c(m-3) - 0.13249937023e11/0.48148892736e11 * c(m-1) + 0.2e1/0.17e2 * c(m) + 0.2083938599e10/0.8024815456e10 * c(m-2) -0.6025413881e10/0.21126554976e11 * c(m-2) - 0.59e2/0.68e2 * c(m) - 0.537416663e9/0.7042184992e10 * c(m-1) 0.3e1/0.3136e4 * c(m-3) + 0.27010400129e11/0.345067064608e12 * c(m-2) + 0.234566387291e12/0.690134129216e12 * c(m-1) + 0.12e2/0.17e2 * c(m);
126 ]; 93 ];
127 94
128 M = M/h; 95 M = M/h;
129 D2 = HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m'); 96 D2 = HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m');
130 end 97 end
131 D2 = @D2_fun; 98 D2 = @D2_fun;
149 M4(m-5:m,m-5:m) = rot90(M4_U, 2); 116 M4(m-5:m,m-5:m) = rot90(M4_U, 2);
150 M4 = 1/h^3*M4; 117 M4 = 1/h^3*M4;
151 118
152 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); 119 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
153 end 120 end
121
122
123 % Alternatives for variable 2nd derivative
124 % % ALTERNATIVES %%%%%%%%%%%%%
125 % % for i=4:m-3
126 % % M(i,i-2:i+2)=[-c(i-1)/0.6e1 + c(i-2)/0.8e1 + c(i)/0.8e1 -c(i-2)/0.6e1 - c(i+1)/0.6e1 - c(i-1)/0.2e1 - c(i)/0.2e1 c(i-2)/0.24e2 + 0.5e1/0.6e1 * c(i-1) + 0.5e1/0.6e1 * c(i+1) + c(i+2)/0.24e2 + 0.3e1/0.4e1 * c(i) -c(i-1)/0.6e1 - c(i+2)/0.6e1 - c(i)/0.2e1 - c(i+1)/0.2e1 -c(i+1)/0.6e1 + c(i)/0.8e1 + c(i+2)/0.8e1;];
127 % % end
128 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
129 % % for i=4:m-3
130 % % M(i,i-2:i+2)= [
131 % % c(i-2)/0.8e1 - c(i-1)/0.6e1 + c(i) /0.8e1 ,
132 % % -c(i-2)/0.6e1 - c(i-1)/0.2e1 - c(i) /0.2e1 - c(i+1)/0.6e1 ,
133 % % c(i-2)/2.4e1 + c(i-1)/1.2e0 + c(i) * 0.3/0.4 + c(i+1)/1.2e0 + c(i+2)/2.4e1 ,
134 % % -c(i-1)/0.6e1 - c(i) /0.2e1 - c(i+1)/0.2e1 - c(i+2)/0.6e1 ,
135 % % c(i) /0.8e1 - c(i+1)/0.6e1 + c(i+2)/0.8e1 ,
136 % % ];
137 % % end
138 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139 % Mm2 = c(r-2)/0.8e1 - c(r-1)/0.6e1 + c(r) /0.8e1 ;
140 % Mm1 = -c(r-2)/0.6e1 - c(r-1)/0.2e1 - c(r) /0.2e1 - c(r+1)/0.6e1 ;
141 % M0 = c(r-2)/2.4e1 + c(r-1)/1.2e0 + c(r) * 0.3/0.4 + c(r+1)/1.2e0 + c(r+2)/2.4e1;
142 % Mp1 = -c(r-1)/0.6e1 - c(r) /0.2e1 - c(r+1)/0.2e1 - c(r+2)/0.6e1;
143 % Mp2 = c(r) /0.8e1 - c(r+1)/0.6e1 + c(r+2)/0.8e1;
144
145 % M(r,:) = spdiags([Mm2 Mm1 M0 Mp1 Mp2],0:2*scheme_radius,length(r),m);
146 % % M(r,:) = spdiags([Mm2 Mm1 M0 Mp1 Mp2],(-2:2)+scheme_radius,M(r,:)); % This is slower
147 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
148 % % %% Somthing is wrong here!!
149 % % Mm2 = c(r-2)/0.8e1 - c(r-1)/0.6e1 + c(r) /0.8e1 ;
150 % % Mm1 = -c(r-2)/0.6e1 - c(r-1)/0.2e1 - c(r) /0.2e1 - c(r+1)/0.6e1 ;
151 % % M0 = c(r-2)/2.4e1 + c(r-1)/1.2e0 + c(r) * 0.3/0.4 + c(r+1)/1.2e0 + c(r+2)/2.4e1;
152 % % Mp1 = -c(r-1)/0.6e1 - c(r) /0.2e1 - c(r+1)/0.2e1 - c(r+2)/0.6e1;
153 % % Mp2 = c(r) /0.8e1 - c(r+1)/0.6e1 + c(r+2)/0.8e1;
154 % % % printSize(M_diag_ind);
155 % % % Mdiags = [Mm2 Mm1 M0 Mp1 Mp2];
156 % % % printSize(Mdiags);
157 % % M(M_diag_ind) = [Mm2 Mm1 M0 Mp1 Mp2]; % This is slightly faster
158 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159 % % Kan man skriva det som en multiplikation av en 3-dim matris?
160 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
161
162
163
164
165
166
167
168
169
170
171