Skip to content

Commit 732cfa0

Browse files
lambdacdmlambdacdm
authored
Add files via upload
1 parent ab4faef commit 732cfa0

File tree

3 files changed

+223
-209
lines changed

3 files changed

+223
-209
lines changed

‎ODE.h‎

Lines changed: 197 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,197 @@
1+
#ifndef ODE_H
2+
#define ODE_H
3+
#include "numerical.h"
4+
using namespace nmr;
5+
6+
namespace ode{
7+
8+
//-----声明部分-----
9+
10+
//常微分方程
11+
template<class DS,class DO> DO DSolve(std::function<DO(DS,DO)>,const pair<DS,DO>&,DS,string);
12+
template <class DS, class DO>DO DSolve(std::function<DO(DS, DO)>, const pair<DS, DO> &, DS, int);
13+
template <class DS, class DO> DO DSolve(std::function<DO(DS, DO)>, const pair<DS, DO> &, DS);
14+
template <class DS, class DO>
15+
std::function<DO(DS)> DSolve(std::function<DO(DS, DO)>, const pair<DS, DO> &, string);
16+
template <class DS, class DO>
17+
std::function<DO(DS)> DSolve(std::function<DO(DS, DO)>, const pair<DS, DO> &);
18+
19+
//常微分方程
20+
template<class DS,class DO,class DP> DO Iteration(const DP &phi,const pair<DS,DO> &initial,DS x,int n)
21+
{
22+
if(n==0)
23+
return get<1>(initial);
24+
DS h= (x - get<0>(initial))/n;//步长
25+
DS xk = get<0>(initial);
26+
DO yk= get<1>(initial);
27+
for (int i = 0; i < n;++i)
28+
{
29+
phi(xk, yk, h);
30+
xk = xk + h;
31+
}
32+
return yk;
33+
}
34+
template<class DS,class DO> DO RKSolve(std::function<DO(DS,DO)> f,const pair<DS,DO>& initial,DS x,int n)
35+
{
36+
auto phi = [=](DS &xk, DO &yk, DS &h) {
37+
DS halfh = h / 2.0;
38+
DS xk_plus_halfh = xk + halfh;
39+
DO k1=f(xk,yk);
40+
DO k2 =f(xk_plus_halfh,yk+halfh*k1);
41+
DO k3 =f(xk_plus_halfh,yk+halfh*k2);
42+
DO k4=f(xk+h,yk+h*k3);
43+
yk = yk + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4);
44+
};
45+
return Iteration(phi, initial, x, n);
46+
}
47+
template<class DS,class DO,class DP,class DF>
48+
DO LinearMultiStep(const DP &phi,const DF &f,const pair<DS,DO> &initial,DS x,int n,int startstep)
49+
{
50+
if(n==0)
51+
return get<1>(initial);
52+
DS h = (x - get<0>(initial))/n;
53+
std::deque<DS> xk={get<0>(initial)};
54+
std::deque<DO> yk={get<1>(initial)};
55+
std::deque<DO> fk={f(xk[0],yk[0])};
56+
for (int i =1; i<=std::min(startstep-1,n);++i)
57+
{
58+
yk.push_back(RKSolve(f,{xk.back(),yk.back()},xk.back()+h,1));
59+
xk.push_back(xk.back()+h);
60+
fk.push_back(f(xk.back(), yk.back()));
61+
}
62+
DO y0adpt = yk.back();
63+
for (int i = startstep; i <=n;++i)
64+
{
65+
phi(xk, yk, fk,h,y0adpt);
66+
xk.push_back(xk.back() + h);
67+
fk.push_back(f(xk.back(), yk.back()));
68+
xk.pop_front();
69+
yk.pop_front();
70+
fk.pop_front();
71+
}
72+
return yk.back();
73+
}
74+
template<class DS,class DO> DO DSolve(std::function<DO(DS,DO)> f,const pair<DS,DO>& initial,DS x,int n,string str)
75+
{
76+
if(str=="Euler"|| str=="forward Euler" || str=="explicit Euler")
77+
{
78+
auto phi = [=](DS &xk, DO &yk, DS &h)
79+
{ yk = yk + h * f(xk, yk);
80+
};
81+
return Iteration(phi, initial, x, n);
82+
}
83+
if(str=="backward Euler" || str=="implicit Euler")
84+
{
85+
auto phi = [=](DS &xk, DO &yk, DS &h) {
86+
DS x_plus_h = xk + h;
87+
DO y=yk+h*f(xk,yk);
88+
for (unsigned i = 0; i <3;++i)
89+
y = yk + h * f(x_plus_h, y);
90+
yk = yk + h * f(x_plus_h, y);
91+
};
92+
return Iteration(phi, initial, x, n);
93+
}
94+
if(str=="trapz" || str=="trapezoidal")
95+
{
96+
auto phi = [=](DS &xk, DO &yk, DS &h) {
97+
DO fk=f(xk,yk);
98+
DS halfh = h / 2.0;
99+
DS x_plus_h = xk + h;
100+
DO y0 = yk + h * fk;
101+
DO y1=yk+halfh*(fk+f(x_plus_h,y0));
102+
yk = yk + halfh * (fk + f(x_plus_h, y1));
103+
};
104+
return Iteration(phi, initial, x, n);
105+
}
106+
if(str=="Heun" || str=="improved Euler")
107+
{
108+
auto phi = [=](DS &xk, DO &yk, DS &h) {
109+
DO k1=f(xk,yk);
110+
DO k2=f(xk+h,yk+h*k1);
111+
yk=yk+h/2.0*(k1+k2);
112+
};
113+
return Iteration(phi, initial, x, n);
114+
}
115+
if(str=="midpoint")
116+
{
117+
auto phi = [=](DS &xk, DO &yk, DS &h) {
118+
DS halfh = h / 2.0;
119+
DO k1=f(xk,yk);
120+
DO k2=f(xk+halfh,yk+halfh*k1);
121+
yk=yk+h*k2;
122+
};
123+
return Iteration(phi, initial, x, n);
124+
}
125+
if(str=="RK4" || str=="Runge-Kutta" || str=="RK")
126+
{
127+
return RKSolve(f, initial, x, n);
128+
}
129+
if(str=="Adams-Bashforth")
130+
{
131+
auto phi = [=](deque<DS> &xk, deque<DO> &yk, deque<DO> &fk, DS &h, DO &y0adpt) {
132+
yk.push_back(yk[1] + 3.0 / 2.0 * h * fk[1] - 1.0 / 2.0 * h * fk[0]);
133+
};
134+
return LinearMultiStep(phi, f, initial, x, n, 2);
135+
}
136+
if(str=="Adams")
137+
{
138+
auto phi = [=](deque<DS> &xk, deque<DO> &yk, deque<DO> &fk, DS &h,DO &y0adpt) {
139+
yk.push_back(yk[3]+h/24.0*(55.0*fk[3]-59.0*fk[2]+37.0*fk[1]-9.0*fk[0]));
140+
};
141+
return LinearMultiStep(phi, f, initial, x, n, 4);
142+
}
143+
if(str=="Milne")
144+
{
145+
auto phi = [=](deque<DS> &xk, deque<DO> &yk, deque<DO> &fk, DS &h,DO &y0adpt) {
146+
yk.push_back(yk[0]+4.0*h/3.0*(2.0*fk[3]-fk[2]+2.0*fk[1]));
147+
};
148+
return LinearMultiStep(phi, f, initial, x, n, 4);
149+
}
150+
if(str=="adaptive Adams")
151+
{
152+
auto phi = [=](deque<DS> &xk, deque<DO> &yk, deque<DO> &fk, DS &h,DO &y0adpt) {
153+
DO y0=yk[3]+h/24.0*(55.0*fk[3]-59.0*fk[2]+37.0*fk[1]-9.0*fk[0]);
154+
DO y1 = y0 + 251.0 * (yk[3] - y0adpt)/270.0;
155+
DO y2 = yk[3] + h / 24.0 * (9.0 * f(xk[3] + h, y1) + 19.0 * fk[3] - 5.0* fk[2] + fk[1]);
156+
yk.push_back(y2-19.0*(y2-y0)/270.0);
157+
y0adpt = y0;
158+
};
159+
return LinearMultiStep(phi, f, initial, x, n, 4);
160+
}
161+
if(str=="adaptive Milne-Hamming")
162+
{
163+
auto phi = [=](deque<DS> &xk, deque<DO> &yk, deque<DO> &fk, DS &h,DO &y0adpt) {
164+
DO y0=yk[0]+4.0*h/3.0*(2.0*fk[3]-fk[2]+2.0*fk[1]);
165+
DO y1=y0 + 112.0*(yk[3]-y0adpt)/121.0;
166+
DO y2 = (9.0 * yk[3] - yk[1]) / 8.0 + 3.0 * h / 8.0 * (f(xk[3] + h, y1)+ 2.0 * fk[3] - fk[2]);
167+
yk.push_back(y2 - 9.0 * (y2 - y0) / 121.0);
168+
y0adpt = y0;
169+
};
170+
return LinearMultiStep(phi, f, initial, x, n, 4);
171+
}
172+
cerr<<"错误：没有定义该方法"<<'\n';
173+
return x;
174+
}
175+
template<class DS,class DO> DO DSolve(std::function<DO(DS,DO)> f,const pair<DS,DO>& initial,DS x,string str)
176+
{
177+
return DSolve(f, initial, x, 100, str);
178+
}
179+
template<class DS,class DO> DO DSolve(std::function<DO(DS,DO)> f,const pair<DS,DO>& initial,DS x,int n)
180+
{
181+
return DSolve(f, initial, x, n,"adaptive Milne-Hamming");
182+
}
183+
template<class DS,class DO> DO DSolve(std::function<DO(DS,DO)> f,const pair<DS,DO>& initial,DS x)
184+
{
185+
return DSolve(f, initial, x, 100);
186+
}
187+
template<class DS,class DO> std::function<DO(DS)> DSolve(std::function<DO(DS,DO)> f,const pair<DS,DO> &initial,string str)
188+
{
189+
return [=](DS x) { return DSolve(f, initial, x, str); };
190+
}
191+
template<class DS,class DO> std::function<DO(DS)> DSolve(std::function<DO(DS,DO)> f,const pair<DS,DO> &initial)
192+
{
193+
return [=](DS x) { return DSolve(f, initial, x); };
194+
}
195+
196+
}
197+
#endif

‎linearalgebra.h‎

Lines changed: 0 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -117,28 +117,6 @@ class Matrix
117117
template <class DB> friend DB Norm(const Matrix<DB>&, string);
118118
template <class DB> friend double Cond(const Matrix<DB>&, double);
119119
template <class DB> friend double Cond(const Matrix<DB>&, string);
120-
//插值
121-
/*template <class DB> friend vector<std::function<DB(DB)>> LagrangeBasis(const Matrix<DB> &);
122-
template <class DB> friend DB LBasisToInterp(const vector<std::function<DB(DB)>>&,const Matrix<DB> &,DB);
123-
template <class DB> friend Matrix<DB> DifferenceQuotientTable(const Matrix<DB> &);
124-
template <class DB> friend Matrix<DB> DividedDifferenceTable(const Matrix<DB> &);
125-
template <class DB> friend DB DifferenceQuotient(const Matrix<DB> &);
126-
template <class DB> friend DB DifferenceQuotient(std::function<DB(DB)>, const Matrix<DB> &);
127-
template <class DB> friend DB DQTableToInterp(const Matrix<DB> &, const Matrix<DB> &,DB);
128-
template <class DB> friend DB DDTableToInterp(const Matrix<DB> &, const Matrix<DB> &,DB);
129-
template <class DB> friend DB _2P3DHermiteInterp(const Matrix<DB> &,int i,int j,DB);
130-
template <class DB> friend DB _3P3DHermiteInterp(const Matrix<DB> &,DB);
131-
template <class DE> friend Matrix<DE> SplineSlope(const Matrix<DE> &);
132-
template <class DE> friend DE Interpolation(const Matrix<DE> &,DE,string);
133-
template <class DE> friend DE Interpolation(const Matrix<DE> &,DE);
134-
template <class DE> friend Matrix<DE> Interpolation(const Matrix<DE> &A, const Matrix<DE> &x, string);
135-
template <class DE> friend std::function<DE(DE)> Interpolation(const Matrix<DE> &,string);
136-
//拟合与逼近
137-
template <class DB> friend Matrix<DB> PseudoInverse(const Matrix<DB>&);
138-
template <class DB> friend Matrix<DB> LeastSquares(const Matrix<DB>&, const Matrix<DB>&);
139-
template <class DB> friend Polynomial<DB> PolyFit(const Matrix<DB> &, const vector<Polynomial<DB>> &);
140-
template <class DB> friend Polynomial<DB> PolyFit(const Matrix<DB>&,int);
141-
template <class DB> friend Polynomial<DB> PolyFit(std::function<DB(DB)>,const vector<Polynomial<DB>>&,DB,DB);*/
142120
private:
143121
vector<vector<DM>> value;
144122
int row_num;

0 commit comments

Comments
 (0)