-
Notifications
You must be signed in to change notification settings - Fork 0
/
P2-(e)-BTCS.cpp
81 lines (60 loc) · 2.02 KB
/
P2-(e)-BTCS.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#define EIGEN_STACK_ALLOCATION_LIMIT 0
#include<bits/stdc++.h>
#include<math.h>
#include<Eigen/Eigen>
#include<iomanip>
#define Pi acos(-1) //定義Pi
using namespace std;
using namespace Eigen;
void NormSolErr_BTCS(){
const int ReDeltaX=100; //DeltaX的倒數 變更格子大小請在此更動
double DeltaX=(double)1/ReDeltaX;
const int ReDeltaT=100000; //DeltaT的倒數 變更格子大小請在此更動
double DeltaT=(double)1/ReDeltaT;
double MaxNorm=0;
Matrix<double,ReDeltaX-1,1>NumSol_Old,NumSol_New,ExSol,SolErr,Boundary_Condition;
NumSol_Old.setZero();
NumSol_New.setZero();
ExSol.setZero();
SolErr.setZero();
Boundary_Condition.setZero();
Boundary_Condition(0,0)=1;
Boundary_Condition(ReDeltaX-2,0)=2;
Matrix<double,ReDeltaX-1,ReDeltaX-1> CoeMatrix;
CoeMatrix.setZero();
for(int i=0;i<=ReDeltaX-2;i++){
for(int j=0;j<=ReDeltaX-2;j++){
if(i==j){
CoeMatrix(i,j)=1+2*(DeltaT/pow(DeltaX,2));
}
else if(i==j+1 || i==j-1){
CoeMatrix(i,j)=-1*(DeltaT/pow(DeltaX,2));
}
}
}
for(int i=0;i<=ReDeltaX-2;i++){
NumSol_Old(i,0)=sin(Pi*((i+1)*DeltaX))+((i+1)*DeltaX)+1;
}
double normSolErr=0;
for(int j=1;j<=ReDeltaT;j++){
NumSol_New=CoeMatrix.colPivHouseholderQr().solve(NumSol_Old+(Boundary_Condition*(DeltaT/pow(DeltaX,2))));
for(int i=0;i<=ReDeltaX-2;i++){
ExSol(i,0)=expl(-1*pow(Pi,2)*j*DeltaT)*sin(Pi*((i+1)*DeltaX))+((i+1)*DeltaX)+1;
}
for(int i=0;i<=ReDeltaX-2;i++){
SolErr(i,0)=abs(ExSol(i,0)-NumSol_New(i,0));
if(SolErr(i,0)>MaxNorm)
{
MaxNorm=SolErr(i,0);
}
normSolErr+=pow(SolErr(i,0),2);
}
NumSol_Old=NumSol_New;
}
normSolErr=sqrt(normSolErr/((ReDeltaX+1)*(ReDeltaT+1)));
cout<<"2 Norm: "<<normSolErr<<endl;
cout<<"Max Norm: "<<MaxNorm<<endl;
}
int main(){
NormSolErr_BTCS();
}