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 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
| #include <cassert>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <map>
#include "ufunction.hpp"
#include "Mesh2d.hpp"
#include "RNM.hpp"
#include "GC.hpp"
#include "gmres.hpp"
#include "cputime.h"
#include "SparseMatMap.hpp"
typedef Mesh2 Mesh;
typedef Mesh::Element Element;
typedef Mesh::Rd Rd;
//typedef typename Mesh::Element Element;
//typedef typename Mesh::Rd Rd;
const int nve = Element::nv;
typedef double R;
class Black_Scholes2D_SparseMatrix {
private:
/*** Attributs privés ***/
int K;
R sigma1, sigma2, q, T, r, dt, dK, dq, L;
Mesh Th;
SparseMatrix<R> A, M;
KN<R> x;
KN<R> x_dK;
/*** Fonction qui lit les données du fichier Java **/
void lireDonneesJava() {
ifstream f("sortie.txt");
if (f) {
f >> K >> r >> sigma1 >> sigma2 >> q >> T >> dt >> L;
f.close();
}
else
cerr << "*** Impossible d'ouvrir le fichier ***" << endl;
}
/*** Definition de la fonction u0 ***/
R u0(const R2 & P){return max(K-exp(P.x)-exp(P.y), 0.);}
/*** Definition des conditions limites ***/
R g(const R2 & P){return 0.;}// initialization
/**** fonctinon f ***/
R f(const R2& P){return u0(P)*exp(-r*T);}
/*** Methode qui calcule la contribution sur chaque element ***/
void MatElement(const Element & K, R alpha, R beta, R matK[Element::nv][Element::nv])
{
const int nve = Element::nv;
const int d = Rd::d;
double clilj = 1./((d+2)*(d+1)); // formule magique
Rd Gradw[nve];
K.Gradlambda(Gradw);
R mu1 = pow(sigma1,2)/2 - r;
R mu2 = pow(sigma2,2)/2 - r;
R2 mu(mu1, mu2);
R cgg=K.mesure()*beta, cuu=K.mesure()*alpha*clilj;
for(int i=0; i<nve; ++i)
for(int j=0; j<=i; ++j) {
R k1 = pow(sigma1, 2)*Gradw[i].x + q*sigma1*sigma2*Gradw[i].y;
R k2 = q*sigma1*sigma2*Gradw[i].x + pow(sigma2, 2)*Gradw[i].y;
R2 kappaGradw(k1, k2);
//R aire2 = pow(K.mesure(), 2);
matK[j][i] = matK[i][j] = cgg*((1./2)*(kappaGradw,Gradw[j])+ ((i==j)/(d+1))*(mu, Gradw[j])) + cuu*(1. + (i==j));
/*matK[j][i] = matK[i][j] = cgg*((1./2)*(kappaGradw,Gradw[j])+ ((i==j)*1./((d+1)))*(mu, Gradw[j])) + cuu*(1. + (i==j));*/
}
}
/*** Assemblage de la matrice ***/
template<class MatriceCreuse>
void BuildMatMap(MatriceCreuse &M, const R alpha,const R beta)
{ typedef typename Mesh::Element Element;
typedef typename Mesh::Rd Rd;
int iK[nve];
R matK[nve][nve];
for (int k=0; k<Th.nt; k++){
const Element& K(Th[k]);
for(int i=0; i<nve; i++)
iK[i]=Th(K[i]);
MatElement(K, alpha, beta, matK);
// assemblage
for (int i=0; i<nve; ++i)
for (int j=0; j<nve; ++j)
if(fabs(matK[i][j])>1e-30)
M[make_pair(iK[i],iK[j])] += matK[i][j];
}
}
/*** Méthode qui imprime un vecteur dans le fichier de nom passé en paramètre */
void printX_InFile(ofstream& fileName, KN<R> X) {
for(int it=0; it<Th.nt; it++)
fileName << (R2) Th[it][0] << " " << X[Th(it,0)]<< endl
<< (R2) Th[it][1] << " " << X[Th(it,1)] << endl
<< (R2) Th[it][2] << " " << X[Th(it,2)] << endl
<< (R2) Th[it][0] << " " << X[Th(it,0)] << endl << endl << endl;
}
public:
Black_Scholes2D_SparseMatrix(const char* mesh);
/*** Méthode qui résoud le problème ***/
template<class R,class Mesh,class Matrix>
void solve() {
typedef typename Mesh::Element Element;
typedef typename Mesh::Rd Rd;
const int nve = Element::nv;
A.clear(); // suprime tous les termes de la matrix
M.clear();
BuildMatMap(A, r + 1./dt, 1); // Matrice de rigidité
BuildMatMap(M, 1./dt, 0.); // Matrice de Masse
cout << " Nombre de coef non nulles de A " << A.size() << "\n";
cout << " Nombre de coef non nulles de M " << M.size() << "\n";
KN<R> u0h(Th.nv), u_old(Th.nv);
for (int k=0;k<Th.nt;k++)
{
const Element& K(Th[k]);
for (int i=0; i<nve; i++) {
u0h[Th(k,i)] = u0(K[i]);
}
}
KN<R> b(Th.nv);
KN<R> e(Th.nv);
b = M*u0h;
KN<R> Pii(Th.nv);
// Applciation des ocnditions aux bord grace a la methode TGV(Tres Grande Valeur)
int ii=0;
double tgv =1e30;
for (int i=0; i<Th.nv; i++)
if(Th(i).lab==2 || Th(i).lab==3) //onGamma())
{
ii++;
A[make_pair(i,i)]=tgv;
b[i]=tgv*g(Th(i));
//b_dK[i]=tgv*g(Th(i));
Pii[i]=1./tgv;
}
else
Pii[i]=1./A[make_pair(i,i)];
MatrixDiag<R> P(Pii);
for(int k=0; k<Th.nv; k++)
e[k]=g(Th(k));
cout << " nb sommet on gamma2 and gamma3 = " << ii << endl;
ofstream file_u0("u0.sol");
ofstream file_ufinal("ufinal.sol");
ofstream file_u0_format_FreeFem("u0C++_format_FreeFem.sol");
ofstream file_ufinal_format_FreeFem("ufinalC++_format_FreeFem.sol");
/*** Ecriture de la solution initale dans les fichiers en format FreeFem et Gnuplot***/
printX_InFile(file_u0, u0h);
file_u0_format_FreeFem << u0h;
R cpu0= CPUtime();
R eps;
/*** BOUCLE EN TEMPS ***/
int isConvergent;
u_old=u0h;
for(int t=1; t*dt<=T; t++) {
b=M*u_old;
int nbkylov = 50, nitermax = Th.nv;
KNM<R> H(nbkylov+1, nbkylov+1);
//isConvergent = GradienConjugue(A,P, b,x,Th.nv,eps=1e-10);
isConvergent = GMRES(A, x, b, P, H, nbkylov, nitermax, eps=1e-10);
if(t%100==0){
cout << "Iteration n°" << t << " => ";
cout << " GMRES: nb iter : " << nitermax << " residu : " << eps << " convergence : " << isConvergent << endl;
}
u_old=x;
}
if (!isConvergent)
cout << " erreur Non convergence residu = " << eps << "\n";
/*** Ecriture de la solution final dans les fichiers en format FreeFem et Gnuplot***/
printX_InFile(file_ufinal, x);
file_ufinal_format_FreeFem << x;
R cpu1=CPUtime()-cpu0;
e -= x;
cout << " -- "<< "SparseMatrix" <<" CPU = " << cpu1 << " s , \t"
<< " err = " << e.linfty() <<endl;
}
};
Black_Scholes2D_SparseMatrix::Black_Scholes2D_SparseMatrix(const char* mesh):Th(mesh), A(Th), M(Th), x(Th.nv) {
lireDonneesJava();
}
int main(int argc, char** argv )
{
assert(argc>1); // verifie que le fichier de maillage est passe en argument
Black_Scholes2D_SparseMatrix BS(argv[1]); // prend le ficchier de maillage en argument
BS.solve(); // resoud le probleme
return 0;
} |
Partager