Skip to content

Commit 85793c8

Browse files
lambdacdmlambdacdm
authored
Add files via upload
1 parent 4005396 commit 85793c8

File tree

1 file changed

+52
-52
lines changed

1 file changed

+52
-52
lines changed

‎numerical.h‎

Lines changed: 52 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -2012,6 +2012,14 @@ template <class DB> double Cond(const Matrix<DB> &a, string str)
20122012
{
20132013
return Norm(a, str) * Norm(Inverse(a), str);
20142014
}
2015+
template<class DF> DF InnerProductC(const Matrix<DF> &A,const Matrix<DF> &B)
2016+
{
2017+
return Get(Transpose(A)*B,0,0);
2018+
}
2019+
template<class DF> DF InnerProductR(const Matrix<DF> &A,const Matrix<DF> &B)
2020+
{
2021+
return Get(A*Transpose(B),0,0);
2022+
}
20152023

20162024
//插值
20172025
template <class DB> vector<std::function<DB(DB)>> LagrangeBasis(const Matrix<DB> &A)
@@ -2043,7 +2051,7 @@ template <class DB> vector<std::function<DB(DB)>> LagrangeBasis(const Matrix<DB>
20432051
template <class DB> DB LBasisToInterp(const vector<std::function<DB(DB)>> &lb,const Matrix<DB> &A,DB x)
20442052
{
20452053
int n = RowSize(A);
2046-
DB s = lb[0](x) * A.value[0][1];
2054+
DB s = lb[0](x) * A(0,1);
20472055
for (int i = 1; i < n;++i)
20482056
{
20492057
s += lb[i](x) * A(i,1);
@@ -2068,7 +2076,7 @@ template <class DB> Matrix<DB> DifferenceQuotientTable(const Matrix<DB> & A)
20682076
template <class DB> Matrix<DB> DividedDifferenceTable(const Matrix<DB> &A)
20692077
{
20702078
int data_num= RowSize(A);
2071-
int der_deg = A.column_num-1;
2079+
int der_deg = ColumnSize(A)-1;
20722080
int n = data_num * der_deg;
20732081
int factorial = 1;
20742082
vector<DB> z(n);
@@ -2094,7 +2102,7 @@ template <class DB> Matrix<DB> DividedDifferenceTable(const Matrix<DB> &A)
20942102
template <class DB> DB DifferenceQuotient(const Matrix<DB> &A)
20952103
{
20962104
int n = RowSize(A);
2097-
return DifferenceQuotientTable(A).value[n-1][n-1];
2105+
return (DifferenceQuotientTable(A))(n-1,n-1);
20982106
}
20992107
template <class DB> DB DifferenceQuotient(std::function<DB(DB)> f, const Matrix<DB> &x)
21002108
{
@@ -2113,22 +2121,22 @@ template <class DB> DB DQTableToInterp(const Matrix<DB> &A, const Matrix<DB> &ta
21132121
DB s = table(n-1,n-1);
21142122
for (int i = n - 2; i >= 0;--i)
21152123
{
2116-
s = table.value[i][i] + (x -A(i,0)) * s;
2124+
s = table(i,i) + (x -A(i,0)) * s;
21172125
}
21182126
return s;
21192127
}
21202128
template <class DB> DB DDTableToInterp(const Matrix<DB> &A, const Matrix<DB> &table,DB x)
21212129
{
21222130
int data_num = RowSize(A);
2123-
int der_deg=A.column_num-1;
2131+
int der_deg=ColumnSize(A)-1;
21242132
int n = data_num * der_deg;
21252133
vector<DB> nest(n);
21262134
nest[0] = DB(1);
21272135
for(int i=0;i<n-1;++i)
2128-
nest[i+1] = nest[i] * (x - A.value[i/der_deg][0]);
2136+
nest[i+1] = nest[i] * (x - A(i/der_deg,0));
21292137
DB s=DB(0);
21302138
for (int i = 0;i<n;++i)
2131-
s += nest[i] * table.value[i][i];
2139+
s += nest[i] * table(i,i);
21322140
return s;
21332141
}
21342142
template <class DE> DE _2P3DHermiteInterp(const Matrix<DE> &M,int i,int j,DE x)
@@ -2137,17 +2145,17 @@ template <class DE> DE _2P3DHermiteInterp(const Matrix<DE> &M,int i,int j,DE x)
21372145
DE l1x = (x - M(i,0)) / (M(j,0) - M(i,0));
21382146
DE l0x2 = l0x * l0x;
21392147
DE l1x2 = l1x * l1x;
2140-
return M(i,1) * (1 + 2 * l1x) * l0x2 + M.value[j][1] * (1 + 2 * l0x) * l1x2 +
2141-
M.value[i][2] * (x - M(i,0)) * l0x2 + M.value[j][2] * (x - M(j,0)) *l1x2;
2148+
return M(i,1) * (1 + 2 * l1x) * l0x2 + M(j,1) * (1 + 2 * l0x) * l1x2 +
2149+
M(i,2) * (x - M(i,0)) * l0x2 + M(j,2) * (x - M(j,0)) *l1x2;
21422150
}
21432151
template <class DE> DE _3P3DHermiteInterp(const Matrix<DE> &M,DE x)
21442152
{
21452153
std::function<DE(DE)> L2 = [=](DE t) {
21462154
return DQTableToInterp(M, DifferenceQuotientTable(M),t);
21472155
};
2148-
DE k = (M.value[1][2]-D(L2,M.value[1][0]))/
2149-
((M.value[1][0]-M.value[0][0])*(M.value[1][0]-M.value[2][0]));
2150-
return L2(x) + k * (x - M.value[0][0]) * (x - M.value[1][0]) * (x - M.value[2][0]);
2156+
DE k = (M(1,2)-D(L2,M(1,0)))/
2157+
((M(1,0)-M(0,0))*(M(1,0)-M(2,0)));
2158+
return L2(x) + k * (x - M(0,0)) * (x - M(1,0)) * (x - M(2,0));
21512159
}
21522160
template <class DE> Matrix<DE> SplineSlope(const Matrix<DE> &A)
21532161
{
@@ -2156,21 +2164,21 @@ template <class DE> Matrix<DE> SplineSlope(const Matrix<DE> &A)
21562164
vector<DE> lambda(n);
21572165
vector<DE> mu(n);
21582166
Matrix<DE> g(n, 1);
2159-
h[0] = A.value[1][0] - A.value[0][0];
2167+
h[0] = A(1,0) - A(0,0);
21602168
lambda[0]=DE(0);
21612169
mu[0] = DE(1);
2162-
g.value[0][0] = 3 * (A.value[1][1]-A.value[0][1]) / h[0];
2170+
g(0,0) = 3 * (A(1,1)-A(0,1)) / h[0];
21632171
for (int k = 1; k<= n - 2;++k)
21642172
{
2165-
h[k]=A.value[k+1][0]-A.value[k][0];
2173+
h[k]=A(k+1,0)-A(k,0);
21662174
lambda[k] = h[k] / (h[k] + h[k - 1]);
21672175
mu[k]=1-lambda[k];
2168-
g.value[k][0]=3*(lambda[k]*(A.value[k][1]-A.value[k-1][1])/h[k-1]
2169-
+mu[k]*(A.value[k+1][1]-A.value[k][1])/h[k]);
2176+
g(k,0)=3*(lambda[k]*(A(k,1)-A(k-1,1))/h[k-1]
2177+
+mu[k]*(A(k+1,1)-A(k,1))/h[k]);
21702178
}
21712179
lambda[n-1]=DE(1);
21722180
mu[n-1]=DE(0);
2173-
g.value[n - 1][0] = 3 * (A.value[n - 1][1] - A.value[n - 2][1]) / h[n - 2];
2181+
g(n-1,0) = 3 * (A(n-1,1) - A(n-2,1)) / h[n - 2];
21742182
Matrix<DE> m=TridiagonalSolve({lambda, vector<DE>(n, DE(2)), mu}, g);
21752183
return m;
21762184
}
@@ -2181,7 +2189,7 @@ template <class DE> DE Interpolation(const Matrix<DE> &A,DE x,string str)
21812189
cerr << "错误：至少需要两个点才能插值" << '\n';
21822190
return x;
21832191
}
2184-
if(A.column_num<2)
2192+
if(ColumnSize(A)<2)
21852193
{
21862194
cerr<<"错误：传入的数据矩阵至少两列"<<'\n';
21872195
return x;
@@ -2208,15 +2216,15 @@ template <class DE> DE Interpolation(const Matrix<DE> &A,DE x,string str)
22082216
int k= 1;
22092217
while(k<n)
22102218
{
2211-
if(x<A.value[k][0])
2219+
if(x<A(k,0))
22122220
{
2213-
return A.value[k - 1][1] * (x - A.value[k][0]) / (A.value[k - 1][0] - A.value[k][0]) +
2214-
A.value[k][1] * (x - A.value[k-1][0]) / (A.value[k][0] - A.value[k - 1][0]);
2221+
return A(k-1,1) * (x - A(k,0)) / (A(k-1,0) - A(k,0)) +
2222+
A(k,1) * (x - A(k-1,0)) / (A(k,0) - A(k-1,0));
22152223
}
22162224
++k;
22172225
}
2218-
return A.value[n-2][1]*(x-A.value[n-2][0])/(A.value[n-2][0]-A.value[n-1][0])+
2219-
A.value[n-1][1]*(x-A.value[n-2][0])/(A.value[n-1][0]-A.value[n-2][0]);
2226+
return A(n-2,1)*(x-A(n-2,0))/(A(n-2,0)-A(n-1,0))+
2227+
A(n-1,1)*(x-A(n-2,0))/(A(n-1,0)-A(n-2,0));
22202228
}
22212229
if(str=="pchip" || str=="piecewise Hermite")
22222230
{
@@ -2238,23 +2246,23 @@ template <class DE> DE Interpolation(const Matrix<DE> &A,DE x,string str)
22382246
Matrix<DE> B(2, 3);
22392247
while(k<n)
22402248
{
2241-
if(x<A.value[k][0])
2249+
if(x<A(k,0))
22422250
{
22432251
for (int i = 0; i < 2;++i)
22442252
{
2245-
B.value[i][2] = m.value[k - 1 + i][0];
2253+
B(i,2) = m(k-1+i,0);
22462254
for (int j = 0; j < 2;++j)
2247-
B(i,j) = A.value[k - 1 + i][j];
2255+
B(i,j) = A(k-1+i,j);
22482256
}
22492257
return _2P3DHermiteInterp(B, 0, 1, x);
22502258
}
22512259
++k;
22522260
}
22532261
for (int i = 0; i < 2;++i)
22542262
{
2255-
B.value[i][2] = m.value[n-2+i][0];
2263+
B(i,2) = m(n-2+i,0);
22562264
for (int j = 0; j < 2;++j)
2257-
B(i,j) = A.value[n-2+ i][j];
2265+
B(i,j) = A(n-2+i,j);
22582266
}
22592267
return _2P3DHermiteInterp(B, 0, 1, x);
22602268
}
@@ -2263,7 +2271,7 @@ template <class DE> DE Interpolation(const Matrix<DE> &A,DE x,string str)
22632271
}
22642272
template <class DE> DE Interpolation(const Matrix<DE> &A,DE x)
22652273
{
2266-
if(A.column_num<=2)
2274+
if(ColumnSize(A)<=2)
22672275
return Interpolation(A,x,"linear");
22682276
return Interpolation(A, x, "piecewise Hermite");
22692277
}
@@ -2274,7 +2282,7 @@ template <class DE> Matrix<DE> Interpolation(const Matrix<DE> &A,const Matrix<DE
22742282
cerr << "错误：至少需要两个点才能插值" << '\n';
22752283
return x;
22762284
}
2277-
if(A.column_num<2)
2285+
if(ColumnSize(A)<2)
22782286
{
22792287
cerr<<"错误：传入的数据矩阵至少两列"<<'\n';
22802288
return x;
@@ -2318,19 +2326,19 @@ template <class DE> Matrix<DE> Interpolation(const Matrix<DE> &A,const Matrix<DE
23182326
int j = 0;
23192327
while(j<RowSize(x) && k<n)
23202328
{
2321-
if(x(j,0)<A.value[k][0])
2329+
if(x(j,0)<A(k,0))
23222330
{
2323-
s(j,0)=A.value[k - 1][1] * (x(j,0)- A.value[k][0]) / (A.value[k - 1][0] - A.value[k][0]) +
2324-
A.value[k][1] * (x(j,0) - A.value[k-1][0]) / (A.value[k][0] - A.value[k - 1][0]);
2331+
s(j,0)=A(k-1,1) * (x(j,0)- A(k,0)) / (A(k-1,0) - A(k,0)) +
2332+
A(k,1) * (x(j,0) - A(k-1,0)) / (A(k,0) - A(k-1,0));
23252333
++j;
23262334
continue;
23272335
}
23282336
++k;
23292337
}
23302338
for (int t = j; t < RowSize(x);++t)
23312339
{
2332-
s.value[t][0] =A.value[n-2][1]*(x.value[t][0]-A.value[n-2][0])/(A.value[n-2][0]-A.value[n-1][0])+
2333-
A.value[n-1][1]*(x.value[t][0]-A.value[n-2][0])/(A.value[n-1][0]-A.value[n-2][0]);
2340+
s(t,0) =A(n-2,1)*(x(t,0)-A(n-2,0))/(A(n-2,0)-A(n-1,0))+
2341+
A(n-1,1)*(x(t,0)-A(n-2,0))/(A(n-1,0)-A(n-2,0));
23342342
}
23352343
return s;
23362344
}
@@ -2351,7 +2359,7 @@ template <class DE> Matrix<DE> Interpolation(const Matrix<DE> &A,const Matrix<DE
23512359
++i;
23522360
}
23532361
for (int t = j; t < RowSize(x);++t)
2354-
s.value[t][0] = _2P3DHermiteInterp(A,n-2,n-1, x.value[t][0]);
2362+
s(t,0) = _2P3DHermiteInterp(A,n-2,n-1, x(t,0));
23552363
return s;
23562364
}
23572365
if(str=="spline")
@@ -2364,13 +2372,13 @@ template <class DE> Matrix<DE> Interpolation(const Matrix<DE> &A,const Matrix<DE
23642372
Matrix<DE> B(2, 3);
23652373
while(j<RowSize(x) && k<n)
23662374
{
2367-
if(x(j,0)<A.value[k][0])
2375+
if(x(j,0)<A(k,0))
23682376
{
23692377
for (int i = 0; i < 2;++i)
23702378
{
2371-
B.value[i][2] = m.value[k - 1 + i][0];
2379+
B(i,2) = m(k-1+i,0);
23722380
for (int j = 0; j < 2;++j)
2373-
B(i,j) = A.value[k - 1 + i][j];
2381+
B(i,j) = A(k-1+i,j);
23742382
}
23752383
s(j,0) = _2P3DHermiteInterp(B, 0, 1, x(j,0));
23762384
++j;
@@ -2382,11 +2390,11 @@ template <class DE> Matrix<DE> Interpolation(const Matrix<DE> &A,const Matrix<DE
23822390
{
23832391
for (int i = 0; i < 2;++i)
23842392
{
2385-
B.value[i][2] = m.value[n-2+i][0];
2393+
B(i,2) = m(n-2+i,0);
23862394
for (int j = 0; j < 2;++j)
2387-
B(i,j) = A.value[n-2+ i][j];
2395+
B(i,j) = A(n-2+i,j);
23882396
}
2389-
s.value[t][0]=_2P3DHermiteInterp(B, 0, 1, x.value[t][0]);
2397+
s(t,0)=_2P3DHermiteInterp(B, 0, 1, x(t,0));
23902398
}
23912399
return s;
23922400
}
@@ -2907,14 +2915,6 @@ template<> double Sqrt<double>(double c)
29072915
}
29082916

29092917
//拟合与逼近
2910-
template<class DF> DF InnerProductC(const Matrix<DF> &A,const Matrix<DF> &B)
2911-
{
2912-
return Get(Transpose(A)*B,0,0);
2913-
}
2914-
template<class DF> DF InnerProductR(const Matrix<DF> &A,const Matrix<DF> &B)
2915-
{
2916-
return Get(A*Transpose(B),0,0);
2917-
}
29182918
template<class DI,class DF,class DG> DI InnerProduct(DF f,DG g,DI a,DI b)
29192919
{
29202920
std::function<DI(DI)> h = [&](DI x) { return f(x) * g(x); };
@@ -3033,7 +3033,7 @@ template <class DB> Polynomial<DB> PolyFit(const Matrix<DB> &xy,int deg)
30333033
A(i,0) = 1;
30343034
for(int j=1;j<=deg;++j)
30353035
{
3036-
A(i,j) = A.value[i][j - 1] * xy(i,0);
3036+
A(i,j) = A(i,j-1) * xy(i,0);
30373037
}
30383038
}
30393039
const Matrix<DB> &c=LeastSquares(A, b);

0 commit comments

Comments
 (0)