Skip to content

Commit 6c0dfb0

Browse files
lambdacdmlambdacdm
committed
adjust folders
1 parent 4a16da6 commit 6c0dfb0

File tree

15 files changed

+6500
-6486
lines changed

15 files changed

+6500
-6486
lines changed

‎mathalgorithm参考文档.pdf‎ renamed to ‎doc/mathalgorithm参考文档.pdf‎

File renamed without changes.

‎FFT.h‎ renamed to ‎lib/FFT.h‎

Lines changed: 298 additions & 298 deletions
Large diffs are not rendered by default.

‎ODE.h‎ renamed to ‎lib/ODE.h‎

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

0 commit comments

Comments
 (0)