Главная
База знаний "Allbest"
Программирование, компьютеры и кибернетика
Построение, численное моделирование и анализ одномерной модели гемодинамики
Построение, численное моделирование и анализ одномерной модели гемодинамики
Особенности моделирования гемодинамики. Одномерная модель течения крови в артериях и ее взаимодействия с подвижными стенками. Численное решение дифференциального уравнения с граничными условиями одномерной модели методами прямых и ортогональной прогонки.
| Рубрика | Программирование, компьютеры и кибернетика |
| Вид | курсовая работа |
| Язык | русский |
| Дата добавления | 24.09.2012 |
| Размер файла | 3,9 M |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
end;
sh=t-0.86*fix (t/0.86);
if pr==1
if sh<0.37 c (1) =c (1) +0.001* (max ( (-679.328*sh^2+502.703*sh+7),70));
else
c (1) =c (1) +0.001* (max ( (387.339*sh^2-666.222*sh+293.476),70));
end; end;
if pr==2
if sh<0.37 c (1) =c (1) +0.01* (max ( (-679.328*sh^2+502.703*sh+7),70));
else
c (1) =c (1) +0.01* (max ( (387.339*sh^2-666.222*sh+293.476),70));
end; end;
if pr==3
if sh<0.37 c (1) =c (1) +0.1* (max ( (-679.328*sh^2+502.703*sh+7),70));
else
c (1) =c (1) +0.1* (max ( (387.339*sh^2-666.222*sh+293.476),70));
end; end;
if pr==4
if sh<0.37 c (1) =c (1) + (max ( (-679.328*sh^2+502.703*sh+7),70));
else
c (1) =c (1) + (max ( (387.339*sh^2-666.222*sh+293.476),70));
end; end;
L= [a (1),b (1),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,a (2),b (2),0,0,-a (4),-b (4),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
…
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,a (9),b (9),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
l= [c (1); c (2) - c (4); c (2) - c (5); 0; c (18) - c (26); c (18) - c (27); 0; c (28) - c (34); c (28) - c (35); 0; c (36) +9; c (37) - c (38); c (37) - c (39); 0; c (40) +9; c (41) - c (42); c (41) - c (43); 0; c (44) - c (46); c (44) - c (47); 0; c (48) +9; c (49) +9; c (45) +9; c (50) - c (52); c (50) - c (53); 0; c (54) +9; c (55) +9; c (51) +9; c (29) - c (30); c (29) - c (31); 0; c (32) +9; c (33) +9; c (19) - c (20); c (19) - c (21); 0; c (22) +9; c (23) - c (24); c (23) - c (25); 0; c (16) +9; c (17) +9; c (3) - c (6); c (3) - c (7); 0; c (10) +9; c (11) - c (12); c (11) - c (13); 0; c (14) +9; c (15) +9; c (8) +9; c (9) +9];
[n_l,n] =size (L);
%*************right boundary condition**************************
for j=1: Nvet,
An (j) =V (1,N+1,j);
Qn (j) =V (2,N+1,j);
a (j) =gamma (j) / (A0 (j) * (sqrt (An (j)) +sqrt (A0 (j))));
b (j) = (rho*Qn (j)) / (2*An (j) ^2);
c (j) =gamma (j) / (sqrt (An (j)) +sqrt (A0 (j)));
end;
R= [a (1),b (1),-a (2),-b (2),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
a (1),b (1),0,0,-a (3),-b (3),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
…
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,-1,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
r= [c (1) - c (2); c (1) - c (3); 0; c (4) - c (18); c (4) - c (19); 0; c (26) +9; c (27) - c (28); c (27) - c (29); 0; c (34) +9; c (35) - c (36); c (35) - c (37); 0; c (38) +9; c (39) - c (40); c (39) - c (41); 0; c (42) - c (44); c (42) - c (45); 0; c (46) - c (48); c (46) - c (49); 0; c (47) +9; c (43) - c (50); c (43) - c (51); 0; c (52) - c (54); c (52) - c (55); 0; c (53) +9; c (30) - c (32); c (30) - c (33); 0; c (31) +9; c (20) +9; c (21) - c (22); c (21) - c (23); 0; c (24) +9; c (25) +9; c (5) - c (16); c (5) - c (17); 0; c (6) - c (10); c (6) - c (11); 0; c (12) +9; c (13) - c (14); c (13) - c (15); 0; c (7) - c (8); c (7) - c (9); 0];
%*************assignment ************************************
n_od=n-n_l; % number of linealy independent solutions homogeneous left boundary condition for one vessel
Z_0=zeros ([n,n_od,N+1]);
Z_f=zeros ([n,N+1]);
Y_0=zeros ([n,n_od,N+1]);
Y_f=zeros ([n,N+1]);
Ort=zeros ([n_od+1,n_od+1,N]); %orthogonalising matrices
U=zeros ([Nfun,N+1,Nvet]); %solution
beta=zeros ([n_od+1,N+1]); %coefficients of return marching
%*************processing of left boundary condition****************
[Q,R_left] =qr (L');
Q_trans=Q';
RR= (R_left (1: n_l,1: n_l)) ';
Y_f (1: n_l,1) =inv (RR) *l;
Y_f (:,1) =Q*Y_f (:,1);
Y_0 (:,:,1) =Q (:,n_l+1: n);
Z_0 (:,:,1) =Y_0 (:,:,1);
Z_f (:,1) =Y_f (:,1);
%**************straight marching********************************
for j=1: Nvet;
q=V (2,1,j);
a=V (1,1,j);
c= (gamma (j) *sqrt (a)) / (2*rho*A0 (j));
G1 (:,:,j) =-znak (j) /tau* [2*q/ (q^2/a-c*a), 1/ (- (q/a) ^2+c);
1,0];
f1 (:,j) =znak (j) * [q* (1/tau+Kr/a) / ( (q/a) ^2-c);
a/tau];
end;
%**************dimensional segmentation loop*********************
for i=1: N,
% vessels loop
for j=1: Nvet
q=V (2, i+1,j);
a=V (1, i+1,j);
c= (gamma (j) *sqrt (a)) / (2*rho*A0 (j));
G2 (:,:,j) =-znak (j) /tau* [2*q/ (q^2/a-c*a), 1/ (- (q/a) ^2+c);
1,0];
f2 (:,j) =znak (j) * [q* (1/tau+Kr/a) / ( (q/a) ^2-c);
a/tau];
% solution of Cauchy problem
del=h (j) /M;
del_mat= (G2 (:,:,j) - G1 (:,:,j)) /M;
del_fun= (f2 (:,j) - f1 (:,j)) /M;
z_0=Z_0 ( (j-1) *Nfun+1: j*Nfun,:, i);
z_f=Z_f ( (j-1) *Nfun+1: j*Nfun, i);
% integration
for jj=1: M,
A=G1 (:,:,j) + (jj-1/2) *del_mat; % intermediate matrix
ff=f1 (:,j) + (jj-1/2) *del_fun; % intermediate right part
del_0=A*z_0*del;
z_0=z_0+del_0; %solution of homogeneous equations
del_f= (A*z_f+ff) *del;
z_f=z_f+del_f; %solution of nonhomogeneous equations
end;
% new value
Y_0 ( (j-1) *Nfun+1: j*Nfun,:, i+1) =z_0;
Y_f ( (j-1) *Nfun+1: j*Nfun, i+1) =z_f;
end;
% orthogonalizaiton
Y= [Y_0 (:,:, i+1),Y_f (:, i+1)];
[Qpr,Rpr] =qr (Y);
Z_f (:, i+1) =Qpr (:,n_od+1) *Rpr (n_od+1,n_od+1);
Rpr (n_od+1,n_od+1) =1;
Z_0 (:,:, i+1) =Qpr (:,1: n_od);
Ort (:,:, i) =Rpr (1: (n_od+1),1: (n_od+1));
end;
%***************inversion of right boundary condition**************
alpha=inv (R*Z_0 (:,:,N+1)) * (r-R*Z_f (:,N+1));
%***************return marching********************************
beta (:,N+1) = [alpha; 1];
Z (:,N+1) = [Z_0 (:,:,N+1),Z_f (:,N+1)] *beta (:,N+1);
if N>1,for i=1: N,
m=N+1-i;
beta (:,m) =inv (Ort (:,:,m)) *beta (:,m+1);
Z (:,m) = [Z_0 (:,:,m),Z_f (:,m)] *beta (:,m);
end;
else
beta (:,1) =inv (Ort (:,:,1)) *beta (:,2);
end;
Z (:,1) = [Y_0 (:,:,1),Y_f (:,1)] *beta (:,1);
%***************vessels loop***********************************
for j=1: Nvet,
U (:,:,j) =Z ( (j-1) *Nfun+1: j*Nfun,:); %solution
end;
%***************reserved operators for output result****************
UUL= [U (:,1,1)];
for i=1: Nvet-2,UUL= [UUL (1: 2*i); U (:,1, i+1)];
end;
UUL= [UUL (1: 2*Nvet-2); U (:,1,Nvet)];
UUR= [U (:,N+1,1)];
for i=1: Nvet-2,UUR= [UUR (1: 2*i); U (:,N+1, i+1)];
end;
UUR= [UUR (1: 2*Nvet-2); U (:,N+1,Nvet)];
WL=L*UUL-l;
WR=R*UUR-r;
Lbound (k) =norm (WL);
Rbound (k) =norm (WR);
%***************move to next step******************************
V=U;
% saving results
% for j=1: Nvet
% OA (k,:,j) =U (1,:,j);
% OQ (k,:,j) =U (2,:,j);
%end;
p (k) =gamma (5) * ( (sqrt (U (1,round (N/2),5)) - sqrt (A0 (5))) /A0 (5));
end;
%**********graphing******************************************
for i=1: Kmax
mas (i) =tau*i;
end;
%F=figure;
%F1=figure;
%F2=figure;
F3=figure;
% graphing middle of 2 vessel area and flux
%figure (F); hold;
%plot (mas,p (:,2),'red'); hold; pause;
%figure (F1); hold;
%plot (mas,p (:,3),'blue'); hold; pause;
% graphing middle of 5 vessel area and flux
%figure (F2); hold;
%plot (mas,p (:,4),'red'); hold; pause;
figure (F3); % hold;
plot (mas,p (:),'blue'); %hold; pause;
end;
disp (max (Lbound));
disp (max (Rbound)); pause;
Energy=1;
2. Пример сгенерированного Java кода
public class HemodynamicsModelSolver
implements TreeTraversal
{
private int [] sign;
protected ArterialBinaryTreeModel atm;
private double [] vesselAreas;
private Matrix [] flux;
public double [] vesselBeta;
private double [] vesselLegths; // static
double [] An, Qn, a, b, c;
Matrix [] OA, OQ;
final int NumberOfUnknownFunction = 2;
int SegmentationForIntegrating = 10, NumberOfSegmentation = 20, rho = 1;
double nu = 0.035, gamma = 265868.077635827404094
double tDelta = 0.02, tMin = 0, tMax = 40;
public void solve (ArterialBinaryTreeModel atm, Matrix [] flux, double [] heartPressure)
{
this. atm = atm;
int NumberOfVessels = atm. size ();
vesselAreas = new double [atm. size ()];
vesselBeta = new double [atm. size ()];
vesselLegths = new double [atm. size ()];
sign = new int [atm. size ()];
atm. traverseTree (this);
for (int i = 0; i < atm. size (); i++)
{
vesselBeta [i] = vesselBeta [i] * 1000;
}
double [] An = new double [NumberOfVessels];
double [] Qn = new double [NumberOfVessels];
double [] a = new double [NumberOfVessels];
double [] b = new double [NumberOfVessels];
double [] c = new double [NumberOfVessels];
Matrix [] Y_0 = new Matrix [NumberOfSegmentation + 1];
Matrix [] Z_0 = new Matrix [NumberOfSegmentation + 1];
Matrix [] Ort = new Matrix [NumberOfSegmentation];
Matrix [] G1 = new Matrix [NumberOfVessels];
Matrix [] G2 = new Matrix [NumberOfVessels];
Matrix [] U = new Matrix [NumberOfVessels];
Matrix [] U0 = new Matrix [NumberOfVessels];
Matrix [] V = new Matrix [NumberOfVessels];
OA = new Matrix [NumberOfVessels];
OQ = new Matrix [NumberOfVessels];
int i, j, k, jj;
int Kmax = (int) Math. round (tMax / tDelta);
double sh = 0.0;
for (i = 0; i < NumberOfSegmentation + 1; i++)
{
Y_0 [i] = new Matrix (NumberOfVessels * 2, NumberOfVessels);
Z_0 [i] = new Matrix (NumberOfVessels * 2, NumberOfVessels);
if (i < NumberOfSegmentation)
{
Ort [i] = new Matrix (NumberOfVessels + 1, NumberOfVessels + 1);
}
}
for (i = 0; i < NumberOfVessels; i++)
{
U [i] = new Matrix (2, NumberOfSegmentation + 1);
U0 [i] = new Matrix (2, NumberOfSegmentation + 1);
V [i] = new Matrix (2, NumberOfSegmentation + 1);
OA [i] = new Matrix (Kmax + 1, NumberOfSegmentation + 1);
OQ [i] = new Matrix (Kmax + 1, NumberOfSegmentation + 1);
G1 [i] = new Matrix (2,2);
G2 [i] = new Matrix (2,2);
}
Matrix L = new Matrix (NumberOfVessels, NumberOfVessels * 2);
Matrix L1 = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);
Matrix Right = new Matrix (NumberOfVessels, NumberOfVessels * 2);
Matrix Q = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);
Matrix Q1 = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);
Matrix R_left = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);
Matrix Rr = new Matrix (NumberOfVessels, NumberOfVessels);
Matrix Y_f = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);
Matrix Z_f = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);
Matrix z_f = new Matrix (2, 1);
Matrix z_0 = new Matrix (2, NumberOfVessels);
Matrix Y = new Matrix (NumberOfVessels * 2, NumberOfVessels + 1);
Matrix Y1 = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);
Matrix Qpr = new Matrix (NumberOfVessels * 2, NumberOfVessels * 2);
Matrix Rpr = new Matrix (NumberOfVessels + 1, NumberOfVessels + 1);
Matrix f1 = new Matrix (2, NumberOfVessels);
Matrix f2 = new Matrix (2, NumberOfVessels);
Matrix ff = new Matrix (2, 1);
Matrix del_mat = new Matrix (2,2);
Matrix del_fun = new Matrix (2, 1);
Matrix A_ = new Matrix (2,2);
Matrix alpha = new Matrix (NumberOfVessels, 1);
Matrix Beta = new Matrix (NumberOfVessels + 1, NumberOfSegmentation + 1);
Matrix Z = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);
Matrix Z_ = new Matrix (NumberOfVessels * 2, NumberOfVessels + 1);
Matrix r = new Matrix (NumberOfVessels, 1);
Matrix l = new Matrix (NumberOfVessels, 1);
Matrix Zero = new Matrix (NumberOfVessels * 2, NumberOfVessels + 1);
Matrix ZeroY_f = new Matrix (NumberOfVessels * 2, NumberOfSegmentation + 1);
for (i = 0; i < NumberOfSegmentation + 1; i++)
{
for (j = 0; j < NumberOfVessels; j++)
{
V [j]. set (0, i, vesselAreas [j]);
}
}
for (j = 0; j < NumberOfVessels; j++)
{
OA [j]. setMatrix (0, 0, 0, NumberOfSegmentation, V [j]. getMatrix (0, 0, 0, NumberOfSegmentation));
OQ [j]. setMatrix (0, 0, 0, NumberOfSegmentation, V [j]. getMatrix (1, 1, 0, NumberOfSegmentation));
}
for (k = 0; k < Kmax; k++)
{
tMin += tDelta;
// experiment repeat marching
for (int k2 = 0; k2 < 1; k2++)
{
for (j = 0; j < NumberOfVessels; j++)
{
if (k2 == 0)
{
An [j] = V [j]. get (0, 0);
Qn [j] = V [j]. get (1, 0);
}
else
{
An [j] = U0 [j]. get (0, 0);
Qn [j] = U0 [j]. get (1, 0);
}
a [j] = vesselBeta [j] / (vesselAreas [j] * (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j])));
b [j] = (rho * Qn [j]) / (2 * An [j] * An [j]);
c [j] = vesselBeta [j] / (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j]));
}
if (k < Kmax / 2)
{
System. out. println ("tMin " + tMin);
}
c [0] +=heartPressure [k] * (1 - Math. exp ( - (tMin) * (tMin) / 100));
double [] flux1 = new double [NumberOfVessels];
int col = flux [0]. getColumnDimension ();
int fluxLength = flux. length;
if (k < Kmax / 2)
{
System. out. println ("a " + a [0] + " b " + b [0] + " c " + c [0]);
}
// resolve left matrix {7}
MatrixResolver leftResolver = new MatrixResolver (atm, a, b, c, true, vesselAreas, flux1, tMin);
Matrix leftMatrix = leftResolver. resolveMatrix ();
L = leftResolver. getMatrix (leftMatrix);
l = leftResolver. getVector (leftMatrix);
for (j = 0; j < NumberOfVessels; j++)
{
if (k2 == 0)
{
An [j] = V [j]. get (0, NumberOfSegmentation);
Qn [j] = V [j]. get (1, NumberOfSegmentation);
}
else
{
An [j] = U0 [j]. get (0, NumberOfSegmentation);
Qn [j] = U0 [j]. get (1, NumberOfSegmentation);
}
a [j] = vesselBeta [j] / (vesselAreas [j] * (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j])));
b [j] = (rho * Qn [j]) / (2 * An [j] * An [j]);
c [j] = vesselBeta [j] / (Math. sqrt (An [j]) + Math. sqrt (vesselAreas [j]));
}
// resolve right matrix {8}
MatrixResolver rightResolver = new MatrixResolver (atm, a, b, c, false, vesselAreas, flux1, tMin);
Matrix rightMatrix = rightResolver. resolveMatrix ();
Right = rightResolver. getMatrix (rightMatrix);
r = rightResolver. getVector (rightMatrix);
// experiment repeat marching was here
System. out. println ("k2 " + k2);
// null matrixes Y_f,Z_f,beta, Z_0, Y_0, Ort, U
Y_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfSegmentation, ZeroY_f);
Z_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfSegmentation, ZeroY_f);
Beta. setMatrix (0, NumberOfVessels, 0, NumberOfSegmentation, ZeroY_f. getMatrix (0, NumberOfVessels, 0,NumberOfSegmentation));
for (int i1 = 0; i1 < NumberOfSegmentation + 1; i1++)
{
Z_0 [i1]. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Zero. getMatrix (0, NumberOfVessels * 2 - 1, 0,NumberOfVessels - 1));
Y_0 [i1]. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Zero. getMatrix (0, NumberOfVessels * 2 - 1, 0,NumberOfVessels - 1));
if (i1 < NumberOfSegmentation)
{
Ort [i1]. setMatrix (0, NumberOfVessels, 0, NumberOfVessels, Zero
. getMatrix (0, NumberOfVessels, 0, NumberOfVessels));
}
}
for (int i2 = 0; i2 < NumberOfVessels; i2++)
{
U [i2]. setMatrix (0, 1, 0, NumberOfSegmentation, ZeroY_f. getMatrix (0, 1, 0, NumberOfSegmentation));
}
// end of null matrixes
L1. setMatrix (0, 2 * NumberOfVessels - 1, 0, NumberOfVessels - 1, L. transpose ());
Q = new QRDecomposition (L1). getQFull ();
R_left = new QRDecomposition (L. transpose ()). getR ();
Rr = (R_left. getMatrix (0, NumberOfVessels - 1, 0, NumberOfVessels - 1)). transpose ();
Y_f. setMatrix (0, NumberOfVessels - 1, 0, 0, (Rr. inverse ()). times (l));
Y_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, 0, Q. times (Y_f. getMatrix (0, NumberOfVessels * 2 - 1, 0, 0)));
Y_0 [0]. setMatrix (0, 2 * NumberOfVessels - 1, 0, NumberOfVessels - 1, Q. getMatrix (0, NumberOfVessels * 2 - 1,NumberOfVessels, NumberOfVessels * 2 - 1));
Z_0 [0] = Y_0 [0];
Z_f. setMatrix (0, NumberOfVessels * 2 - 1, 0, 0, Y_f. getMatrix (0, NumberOfVessels * 2 - 1, 0, 0));
for (j = 0; j < NumberOfVessels; j++) // straight marching
{
double q, a_, q_k, a_k;
if (k2 == 0)
{
q = V [j]. get (1, 0);
a_ = V [j]. get (0, 0);
q_k = 0;
a_k = 0;
}
else
{
q = U0 [j]. get (1, 0);
a_ = U0 [j]. get (0, 0);
q_k = V [j]. get (1, 0);
a_k = V [j]. get (0, 0);
}
double c_ = (vesselBeta [j] * Math. sqrt (a_)) / (2 * rho * vesselAreas [j]);
G1 [j]. set (0, 0, ( - sign [j] / tDelta) * 2 * q / (q * q / a_ - c_ * a_));
G1 [j]. set (0, 1, ( - sign [j] / tDelta) / ( - (q / a_) * (q / a_) + c_));
G1 [j]. set (1, 0, - sign [j] / tDelta);
G1 [j]. set (1, 1, 0);
if (k2 == 0)
{
f1. set (0, j, sign [j] * (q * (1/tDelta + (8 * Math. PI * nu) / a_) / ( (q / a_) * (q / a_) - c_)));
f1. set (1, j, sign [j] * a_ / tDelta);
}
else
{
f1. set (0, j, sign [j]
* ( (1/tDelta) * (q_k + (8 * Math. PI * nu) * tDelta * q / a_) / ( (q / a_) * (q / a_) - c_)));
f1. set (1, j, sign [j] * a_k / tDelta);
}
}
for (i = 0; i < NumberOfSegmentation; i++) // dimensional
// segmentation loop
{
for (j = 0; j < NumberOfVessels; j++) // vessel loop
{
double q, a_, q_k, a_k;
if (k2 == 0)
{
q = V [j]. get (1, i + 1);
a_ = V [j]. get (0, i + 1);
q_k = 0;
a_k = 0;
}
else
{
q = U0 [j]. get (1, i + 1);
a_ = U0 [j]. get (0, i + 1);
q_k = V [j]. get (1, i + 1);
a_k = V [j]. get (0, i + 1);
}
double c_ = (vesselBeta [j] * Math. sqrt (a_)) / (2 * rho * vesselAreas [j]);
G2 [j]. set (0, 0, ( - sign [j] / tDelta) * 2 * q / (q * q / a_ - c_ * a_));
G2 [j]. set (0, 1, ( - sign [j] / tDelta) / ( - (q / a_) * (q / a_) + c_));
G2 [j]. set (1, 0, - sign [j] / tDelta);
G2 [j]. set (1, 1, 0);
if (k2 == 0)
{
f2. set (0, j, sign [j] * (q * (1/tDelta + (8 * Math. PI * nu) / a_) / ( (q / a_) * (q / a_) - c_)));
f2. set (1, j, sign [j] * a_ / tDelta);
}
else
{
f2. set (0, j,
sign [j]
* ( (1/tDelta) * (q_k + (8 * Math. PI * nu) * tDelta * q / a_) / ( (q / a_)
* (q / a_) - c_)));
f2. set (1, j, sign [j] * a_k / tDelta);
}
// solution of Cauchy problem
double del = vesselLegths [j] / (NumberOfSegmentation * SegmentationForIntegrating);
del_mat. setMatrix (0, 1, 0, 1, (G2 [j]. minus (G1 [j])). times (1.0/SegmentationForIntegrating));
del_fun. set (0, 0, (f2. get (0, j) - f1. get (0, j)) / SegmentationForIntegrating);
del_fun. set (1, 0, (f2. get (1, j) - f1. get (1, j)) / SegmentationForIntegrating);
z_0 = Z_0 [i]. getMatrix (j * 2, j * 2 + 1, 0, NumberOfVessels - 1);
z_f = Z_f. getMatrix (j * 2, j * 2 + 1, i, i);
for (jj = 0; jj < SegmentationForIntegrating; jj++)
{
A_. setMatrix (0, 1, 0, 1, G1 [j]. plus (del_mat. times (jj + 0.5))); // intermediate
// matrix
ff. setMatrix (0, 1, 0, 0, (f1. getMatrix (0, 1, j, j)). plus (del_fun. times (jj + 0.5))); // intermediate
// right
// part
z_0. plusEquals (A_. times (z_0. times (del)));
z_f. plusEquals ( ( (A_. times (z_f)). plus (ff)). times (del));
}
// new value
Y_0 [i + 1]. setMatrix (j * 2, j * 2 + 1, 0, NumberOfVessels - 1, z_0);
Y_f. setMatrix (j * 2, j * 2 + 1, i + 1, i + 1, z_f);
for (int i3 = 0; i3 < NumberOfVessels; i3++)
{
G1 [i3] = G2 [i3]. copy ();
f1 = f2. copy ();
}
}
// orthogonalizaiton
Y. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Y_0 [i + 1]);
Y. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Y_f. getMatrix (0, NumberOfVessels * 2 - 1,i + 1, i + 1));
Y1. setMatrix (0, 2 * NumberOfVessels - 1, 0, NumberOfVessels, Y);
Qpr = new QRDecomposition (Y1). getQFull ();
Rpr = new QRDecomposition (Y). getR (); // k==1 +
Z_f. setMatrix (0, NumberOfVessels * 2 - 1, i + 1, i + 1, Qpr. getMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels,
NumberOfVessels). times (Rpr. get (NumberOfVessels, NumberOfVessels)));
Rpr. set (NumberOfVessels, NumberOfVessels, 1);
Z_0 [i + 1]. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Qpr. getMatrix (0, NumberOfVessels * 2 - 1,0, NumberOfVessels - 1));
Ort [i]. setMatrix (0, NumberOfVessels, 0, NumberOfVessels, Rpr. getMatrix (0, NumberOfVessels, 0, NumberOfVessels));
}
alpha. setMatrix (0, NumberOfVessels - 1, 0, 0, ( (Right. times (Z_0 [NumberOfSegmentation])). inverse ())
. times (r. minus (Right. times (Z_f. getMatrix (0, NumberOfVessels * 2 - 1, NumberOfSegmentation,
NumberOfSegmentation)))));
// return marching
Beta. setMatrix (0, NumberOfVessels - 1, NumberOfSegmentation, NumberOfSegmentation, alpha);
Beta. set (NumberOfVessels, NumberOfSegmentation, 1);
Z_. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Z_0 [NumberOfSegmentation]. getMatrix (0,NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1));
Z_. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Z_f. getMatrix (0, NumberOfVessels * 2 - 1,NumberOfSegmentation, NumberOfSegmentation));
Z. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfSegmentation, NumberOfSegmentation, Z_. times (Beta. getMatrix (0,NumberOfVessels, NumberOfSegmentation, NumberOfSegmentation)));
if (NumberOfSegmentation > 1)
{
for (i = 0; i < NumberOfSegmentation; i++)
{
int m = NumberOfSegmentation - i;
Beta. setMatrix (0, NumberOfVessels, m - 1, m - 1, (Ort [m - 1]. inverse ()). times (Beta. getMatrix (0,NumberOfVessels, m, m)));
Z_. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Z_0 [m - 1]. getMatrix (0,NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1));
Z_. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Z_f. getMatrix (0,NumberOfVessels * 2 - 1, m - 1, m - 1));
Z. setMatrix (0, NumberOfVessels * 2 - 1, m - 1, m - 1, Z_. times (Beta. getMatrix (0, NumberOfVessels, m - 1,m - 1)));
}
}
else
{
Beta
. setMatrix (0, NumberOfVessels, 0, 0, (Ort [0]. inverse ()). times (Beta
. getMatrix (0, NumberOfVessels, 1, 1)));
}
Z_. setMatrix (0, NumberOfVessels * 2 - 1, 0, NumberOfVessels - 1, Y_0 [0]. getMatrix (0, NumberOfVessels * 2 - 1, 0,NumberOfVessels - 1));
Z_. setMatrix (0, NumberOfVessels * 2 - 1, NumberOfVessels, NumberOfVessels, Y_f. getMatrix (0, NumberOfVessels * 2 - 1, 0,0));
Z. setMatrix (0, NumberOfVessels * 2 - 1, 0, 0, Z_. times (Beta. getMatrix (0, NumberOfVessels, 0, 0)));
// vessels loop
for (j = 0; j < NumberOfVessels; j++)
{
U [j]. setMatrix (0, 1, 0, NumberOfSegmentation, Z. getMatrix (j * 2, j * 2 + 1, 0, NumberOfSegmentation)); // solution
}
for (i = 0; i < NumberOfVessels; i++)
{
U0 [i] = U [i]. copy ();
}
if ( (k2 == 0) && (tMin < 2))
{
break;
}
} // end of experiment repeat marching
for (i = 0; i < NumberOfVessels; i++)
{
V [i] = U [i]. copy ();
}
for (j = 0; j < NumberOfVessels; j++)
{
OA [j]. setMatrix (k + 1, k + 1, 0, NumberOfSegmentation, U [j]. getMatrix (0, 0, 0, NumberOfSegmentation));
OQ [j]. setMatrix (k + 1, k + 1, 0, NumberOfSegmentation, U [j]. getMatrix (1, 1, 0, NumberOfSegmentation));
}
} // end of general loop
}
public void setTDelta (double value)
{
tDelta = value;
}
public double getTDelta ()
{
return tDelta;
}
public void setTMin (double value)
{
tMin = value;
}
public void setTMax (double value)
{
tMax = value;
}
public void setNumberOfSegmentation (int value)
{
NumberOfSegmentation = value;
}
public void setSegmentationForIntegrating (int value)
{
SegmentationForIntegrating = value;
}
public Matrix [] getArea ()
{
return OA;
}
public Matrix [] getFlux ()
{
return OQ;
}
public void visit (Vessel v)
{
int i = v. index;
System. out. println (atm. size ()); // for debug
vesselLegths [i] = v. length;
vesselAreas [i] = v. area;
vesselBeta [i] = v. getBeta ();
if (atm. isEvenNode (v))
sign [i] = 1;
else
sign [i] = - 1;
if (v. left! = null)
visit (v. left);
if (v. right! = null)
visit (v. right);
}
public static class MatrixResolver
implements TreeTraversal
{
public Matrix getMatrix (Matrix m)
{
Matrix resultMatrix = new Matrix (atm. size (), atm. size () * 2);
resultMatrix = (m. getMatrix (0, 2 * atm. size () - 1, 0, atm. size () - 1)). transpose ();
return resultMatrix;
}
public Matrix getVector (Matrix v)
{
Matrix resultVector = new Matrix (atm. size (), 1);
resultVector = (v. getMatrix (2 * atm. size (), 2 * atm. size (), 0, atm. size () - 1)). transpose ();
return resultVector;
}
private double [] a;
private double [] b;
private double [] c;
private Matrix result = null;
private int index = 0;
private boolean left = true;
private double [] vesselAreas;
private int k;
private ArterialBinaryTreeModel atm;
private double [] flux;
private double t;
public MatrixResolver (ArterialBinaryTreeModel atm, double [] a, double [] b, double [] c, boolean left, double [] vesselAreas,
double [] flux, double t)
{
this. atm = atm;
this. a = a;
this. b = b;
this. c = c;
this. left = left;
this. vesselAreas = vesselAreas;
this. flux = flux;
this. t = t;
result = new Matrix (2 * atm. size () + 1, atm. size ());
}
public Matrix resolveMatrix ()
{
if (left)
{
double flux0 = flux [0];
double fluxFin = flux [0];
Matrix singlePressureFirstCondition = resolveSinglePressureEndCondition (0, a, b, c, flux0, fluxFin, t);
result. setMatrix (0, 2 * atm. size (), index, index, singlePressureFirstCondition);
index++;
}
atm. traverseTree (this);
return result;
}
public void visit (Vessel v)
{
vesselAreas [v. index] = v. area;
if (atm. isEvenNode (v) == left)
return;
if (v. left! = null && v. right! = null)
{
Matrix singlePressureConditionLeft = resolveSinglePressureCondition (v, v. left. index, a, b, c);
result. setMatrix (0, 2 * atm. size (), index, index, singlePressureConditionLeft);
index++;
Matrix singlePressureConditionRight = resolveSinglePressureCondition (v, v. right. index, a, b, c);
result. setMatrix (0, 2 * atm. size (), index, index, singlePressureConditionRight);
index++;
Matrix singleFluxCondition = resolveSingleFluxCondition (v. index, v. left. index, v. right. index, a, b, c);
result. setMatrix (0, 2 * atm. size (), index, index, singleFluxCondition);
index++;
}
else
{
double flux0 = flux [0];
double fluxFin = flux [v. index];
Matrix singlePressureEndCondition = resolveSinglePressureEndCondition (v. index, a, b, c, flux0, fluxFin, t);
result. setMatrix (0, 2 * atm. size (), index, index, singlePressureEndCondition);
index++;
}
}
private Matrix resolveSinglePressureCondition (Vessel vessel, int childIndex, double [] a, double [] b, double [] c)
{
Matrix result = new Matrix (2 * atm. size () + 1, 1);
int parentIndex = vessel. index;
result. set (2 * parentIndex, 0, a [parentIndex]);
result. set (2 * parentIndex + 1, 0, b [parentIndex]);
result. set (2 * childIndex, 0, - a [childIndex]);
result. set (2 * childIndex + 1, 0, - b [childIndex]);
result. set (2 * atm. size (), 0, c [parentIndex] - c [childIndex]);
return result;
}
private Matrix resolveSingleFluxCondition (int i, int j, int k, double [] a, double [] b, double [] c)
{
Matrix result = new Matrix (2 * atm. size () + 1, 1);
result. set (2 * i + 1, 0, 1);
result. set (2 * j + 1, 0, - 1);
result. set (2 * k + 1, 0, - 1);
return result;
}
private Matrix resolveSinglePressureEndCondition (int vesselIndex, double [] a, double [] b, double [] c, double flux0,double fluxFin, double t)
{
Matrix result = new Matrix (2 * atm. size () + 1, 1);
result. set (2 * vesselIndex, 0, a [vesselIndex]);
result. set (2 * vesselIndex + 1, 0, b [vesselIndex]);
if (vesselIndex == 0)
{
result. set (2 * atm. size (), 0, c [vesselIndex]);
}
else
{
result. set (2 * atm. size (), 0, c [vesselIndex] + 70 * (1 - Math. exp ( - (t * t) / 100)));
}
return result;
}
}
}
Размещено на Allbest.ru
Подобные документы
Разностная схема решения уравнения теплопроводности. Численное решение уравнения теплопроводности в табличном процессоре Microsoft Ехсеl и в пакете математических расчётов MathCAD. Расчёт методом прогонки. Изменение пространственной координаты.
дипломная работа [248,4 K], добавлен 15.03.2014Разработка программы на языке С++ для решения дифференциального уравнения Лапласа в прямоугольной области методом сеток. Численное решение задачи Дирихле для уравнения Лапласа, построение сетки и итерационного процесса. Листинг и результат программы.
курсовая работа [307,5 K], добавлен 30.04.2012Решение задачи Коши для дифференциального уравнения методом Рунге-Кутта и Адамса с автоматическим выбором шага и заданным шагом. Интерполирование табличной функции. Численное решение системы линейных алгебраических уравнений методами простой итерации.
методичка [35,8 K], добавлен 15.03.2009Численное решение задачи Коши для обыкновенного дифференциального уравнения первого и второго порядка методом Эйлера и Рунге-Кутты и краевой задачи для ОДУ второго порядка с применением пакета MathCad, электронной таблицы Excel и программы Visual Basic.
курсовая работа [476,2 K], добавлен 14.02.2016Математическое описание алгоритмов схемы и операций для уравнения Лапласа. Изучение разностной схемы "крест" для нахождения численного решения эллиптического уравнения, задача Дирихле. Использование указателей в среде Matlab для решений методом Гаусса.
дипломная работа [859,3 K], добавлен 23.10.2014Схема электрической цепи (источник переменного тока, катушка индуктивности, конденсатор, набор резисторов и ключ). Вывод системы дифференциальных уравнений. Численное интегрирование (методы левых и средних прямоугольников). Блок-схемы и программные коды.
курсовая работа [1,7 M], добавлен 09.06.2012Разработка прикладного программного обеспечения для решения расчетных задач для компьютера. Численное интегрирование - вычисление значения определённого интеграла. Проектирование алгоритма численного метода. Тестирование работоспособности программы.
курсовая работа [1,1 M], добавлен 03.08.2011Расчет тепловой схемы с применением методов математического моделирования. Разработка алгоритма реализации модели. Составление программы для ПЭВМ, ее отладка и тестирование. Проведение численного исследования и параметрическая оптимизация системы.
курсовая работа [2,8 M], добавлен 01.03.2013Построение концептуальной модели системы и ее формализация. Алгоритмизация модели системы и ее машинная реализация. Построение логической схемы модели. Проверка достоверности модели системы. Получение и интерпретация результатов моделирования системы.
курсовая работа [67,9 K], добавлен 07.12.2009Построение концептуальной модели и метод имитационного моделирования. Определение переменных уравнений математической модели и построение моделирующего алгоритма. Описание возможных улучшений системы и окончательный вариант модели с результатами.
курсовая работа [79,2 K], добавлен 25.06.2011
