From df1cbb2067b89732c6601ca9c00f7d8bcf7fb58a Mon Sep 17 00:00:00 2001 From: Steve Policastro Date: Thu, 22 Jul 2021 13:55:09 -0400 Subject: [PATCH] Working on the solver. Neumann BCs are still not working properly. And, I need to re-implement the Dirichlet BCs so that I'm accounting for the boundary nodes correctly... --- Diffusion2D_Library.sln | 6 ++ DiffusionSimulators_1D.cs | 10 +-- DiffusionSimulators_2D_ConstantD.cs | 28 +++---- DiffusionSimulators_2D_MatrixD.cs | 119 ++++++++++++++++------------ RMatrix.cs | 8 +- TridiagonalMatrix.cs | 6 +- log.txt | 14 ++++ 7 files changed, 114 insertions(+), 77 deletions(-) diff --git a/Diffusion2D_Library.sln b/Diffusion2D_Library.sln index e5fd45d..2ffbaa8 100644 --- a/Diffusion2D_Library.sln +++ b/Diffusion2D_Library.sln @@ -5,6 +5,8 @@ VisualStudioVersion = 16.0.31402.337 MinimumVisualStudioVersion = 10.0.40219.1 Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Diffusion2D_Library", "Diffusion2D_Library.csproj", "{0402B501-B293-47CA-AE2E-7939BD274CAE}" EndProject +Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Test_Diffusion2D_Library", "..\Test_Diffusion2D_Library\Test_Diffusion2D_Library.csproj", "{43A63030-4802-4B22-AD0F-2A42A17C6C09}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Any CPU = Debug|Any CPU @@ -15,6 +17,10 @@ Global {0402B501-B293-47CA-AE2E-7939BD274CAE}.Debug|Any CPU.Build.0 = Debug|Any CPU {0402B501-B293-47CA-AE2E-7939BD274CAE}.Release|Any CPU.ActiveCfg = Release|Any CPU {0402B501-B293-47CA-AE2E-7939BD274CAE}.Release|Any CPU.Build.0 = Release|Any CPU + {43A63030-4802-4B22-AD0F-2A42A17C6C09}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {43A63030-4802-4B22-AD0F-2A42A17C6C09}.Debug|Any CPU.Build.0 = Debug|Any CPU + {43A63030-4802-4B22-AD0F-2A42A17C6C09}.Release|Any CPU.ActiveCfg = Release|Any CPU + {43A63030-4802-4B22-AD0F-2A42A17C6C09}.Release|Any CPU.Build.0 = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/DiffusionSimulators_1D.cs b/DiffusionSimulators_1D.cs index b5fb097..0a870be 100644 --- a/DiffusionSimulators_1D.cs +++ b/DiffusionSimulators_1D.cs @@ -111,7 +111,7 @@ public void OneD_FE(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia RVector b = new(n); // Define the A matrix - double nu = (D * dt) / (Math.Pow(dx, 2)); + double nu = D * dt / Math.Pow(dx, 2); double off_d_val = nu; double diag_val = 1 - (2 * nu); TridiagonalMatrix A = new(n, diag_val, off_d_val, off_d_val); @@ -158,8 +158,8 @@ public void OneD_FE(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia break; case BoundaryConditions.Neumann: xnew = A.Dot(xold) + g(position, t * dt, xold); - xnew[0] = (-2 * nu) * Lbc(t * dt); - xnew[n - 1] = (2 * nu) * Rbc(t * dt); + xnew[0] = -2 * nu * Lbc(t * dt); + xnew[n - 1] = 2 * nu * Rbc(t * dt); break; case BoundaryConditions.Mixed: xnew = A.Dot(xold) + b; @@ -184,7 +184,7 @@ public void OneD_BE(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia RVector xnew = new(n); // Define the A matrix - double nu = (D * dt) / (Math.Pow(dx, 2)); + double nu = D * dt / Math.Pow(dx, 2); double off_d_val = -nu; double diag_val = 1 + (2 * nu); TridiagonalMatrix A = new(n, diag_val, off_d_val, off_d_val); @@ -254,7 +254,7 @@ public void OneD_CN(BoundaryCondition_Del Lbc, BoundaryCondition_Del Rbc, Initia RVector bj1 = new(n); // Setup the A matrix - double nu = (D * dt) / (Math.Pow(dx, 2)); + double nu = D * dt / Math.Pow(dx, 2); double off_d_val, diag_val; off_d_val = -nu / 2; diag_val = 1.0 + nu; diff --git a/DiffusionSimulators_2D_ConstantD.cs b/DiffusionSimulators_2D_ConstantD.cs index cc5ff63..f12eba4 100644 --- a/DiffusionSimulators_2D_ConstantD.cs +++ b/DiffusionSimulators_2D_ConstantD.cs @@ -255,7 +255,7 @@ public void Solve(int n_time_steps, int output_interval) RVector CT = new(ncols); RVector CB = new(ncols); - double nu = (d_coeff * dt) / (2 * Math.Pow(dx, 2)); + double nu = d_coeff * dt / (2 * Math.Pow(dx, 2)); // Define the A matrix for the explicit steps double off_d_val = nu; @@ -329,9 +329,9 @@ public void Solve(int n_time_steps, int output_interval) cUpper2 = C_Im2[nrows - 2, ncols - 2]; } //C3 = BCs[1].BoundaryFunction(t * dt, ncols); - CR = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3 - CR[0] = ((-2 * nu) * cLower1) + ((2 * nu) * cLower2); - CR[ncols - 1] = ((-2 * nu) * cUpper1) + ((2 * nu) * cUpper2); + CR = (-2 * nu * C1) + (2 * nu * C2); // + C3 + CR[0] = (-2 * nu * cLower1) + (2 * nu * cLower2); + CR[ncols - 1] = (-2 * nu * cUpper1) + (2 * nu * cUpper2); break; default: break; @@ -365,9 +365,9 @@ public void Solve(int n_time_steps, int output_interval) cUpper2 = C_Im2[nrows - 2, 1]; } //C3 = BCs[1].BoundaryFunction(t * dt, ncols); - CL = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3 - CL[0] = ((-2 * nu) * cLower1) + ((2 * nu) * cLower2); - CL[ncols - 1] = ((-2 * nu) * cUpper1) + ((2 * nu) * cUpper2); + CL = (-2 * nu * C1) + (2 * nu * C2); // + C3 + CL[0] = (-2 * nu * cLower1) + (2 * nu * cLower2); + CL[ncols - 1] = (-2 * nu * cUpper1) + (2 * nu * cUpper2); break; default: break; @@ -404,7 +404,7 @@ public void Solve(int n_time_steps, int output_interval) } //f12 = (tau / 2.0) * (fn + f0); - f12 = (tau / 2.0) * fn; + f12 = tau / 2.0 * fn; // BCs switch (BCs_Functions[0].TypeBC) @@ -412,7 +412,7 @@ public void Solve(int n_time_steps, int output_interval) case ABoundaryCondition.dirichlet: RVector gn = BCs_Functions[0].BoundaryFunction((t + 1) * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed); RVector g0 = BCs_Functions[0].BoundaryFunction(t * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed); - CT = ((0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0))); ////((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu * + CT = (0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0)); ////((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu * break; case ABoundaryCondition.neumann: RVector C1, C2; //,C3; @@ -427,7 +427,7 @@ public void Solve(int n_time_steps, int output_interval) C2 = C_Im2.GetRowVector(nrows - 2); } //C3 = BCs[1].BoundaryFunction(t * dt, ncols); - CT = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3 + CT = (-2 * nu * C1) + (2 * nu * C2); // + C3 break; default: break; @@ -437,7 +437,7 @@ public void Solve(int n_time_steps, int output_interval) case ABoundaryCondition.dirichlet: RVector gn = BCs_Functions[3].BoundaryFunction((t + 1) * dt, BCs_Functions[3].PositionVaries, BCs_Functions[3].PositionFixed); RVector g0 = BCs_Functions[3].BoundaryFunction(t * dt, BCs_Functions[3].PositionVaries, BCs_Functions[3].PositionFixed); - CB = ((0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0))); //((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu * + CB = (0.5 * B.Dot(gn)) + (0.5 * A.Dot(g0)); //((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu * break; case ABoundaryCondition.neumann: RVector C1, C2; //,C3; @@ -452,7 +452,7 @@ public void Solve(int n_time_steps, int output_interval) C2 = C_Im2.GetRowVector(1); } //C3 = BCs[1].BoundaryFunction(t * dt, ncols); - CB = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3 + CB = (-2 * nu * C1) + (2 * nu * C2); // + C3 break; default: break; @@ -496,7 +496,7 @@ public void Solve(int n_time_steps, int output_interval) C2 = C_Im2.GetColVector(ncols - 2); } //C3 = BCs[1].BoundaryFunction(t * dt, ncols); - CR = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3 + CR = (-2 * nu * C1) + (2 * nu * C2); // + C3 break; default: break; @@ -524,7 +524,7 @@ public void Solve(int n_time_steps, int output_interval) C2 = C_Im2.GetColVector(1); } //C3 = BCs[1].BoundaryFunction(t * dt, ncols); - CL = ((-2 * nu) * C1) + ((2 * nu) * C2); // + C3 + CL = (-2 * nu * C1) + (2 * nu * C2); // + C3 break; default: break; diff --git a/DiffusionSimulators_2D_MatrixD.cs b/DiffusionSimulators_2D_MatrixD.cs index f04df23..5d7568b 100644 --- a/DiffusionSimulators_2D_MatrixD.cs +++ b/DiffusionSimulators_2D_MatrixD.cs @@ -78,124 +78,144 @@ public DiffusionSimulators_2D_MatrixD( // ==================================== // Create the solver method matrices // ==================================== + // Total points minus the left/right or up/down boundaries + int nx_less1 = nx - 2; + // nx_less1 - one point because the total entries are one fewer because these are on the off-diagonals + int nx_less2 = nx_less1 - 1; + + int end_idx1 = nx_less2; + int end_idx2 = nx_less2 - 1; + int start_idx1 = 0; + int start_idx2 = 1; + A = new TridiagonalMatrix[nx]; B_row = new TridiagonalMatrix[nx]; B_col = new TridiagonalMatrix[nx]; double nu0 = dt / (2 * Math.Pow(dx, 2)); - RVector nu = new(nx); - RVector off_d_val_l0 = new(nx - 1); - RVector off_d_val_u0 = new(nx - 1); - RVector main0 = new(nx); - RVector off_d_val_l1 = new(nx - 1); - RVector off_d_val_u1 = new(nx - 1); - RVector main1 = new(nx); - RVector off_d_val_l2 = new(nx - 1); - RVector off_d_val_u2 = new(nx - 1); - RVector main2 = new(nx); + RVector nu = new(nx_less1); + RVector off_d_val_l0 = new(nx_less2); + RVector off_d_val_u0 = new(nx_less2); + RVector main0 = new(nx_less1); + RVector off_d_val_l1 = new(nx_less2); + RVector off_d_val_u1 = new(nx_less2); + RVector main1 = new(nx_less1); + RVector off_d_val_l2 = new(nx_less2); + RVector off_d_val_u2 = new(nx_less2); + RVector main2 = new(nx_less1); double D_minus, D_plus; - for (int i = 0; i < nx; i++) + for (int i = 0; i < nx-1; i++) { RVector D_row = D.GetRowVector(i); + double D0 = D_row[0]; + double D1 = D_row[1]; + double Dend = D_row[nx - 1]; + double Dend1 = D_row[nx - 2]; // Working this for diffusion heterogeneous media //================================================ // A Matrix first //================================================ - D_plus = 0.5 * (D_row[1] + D_row[0]); - main0[0] = 1.0 - ((D_plus * nu0)); + #region A + D_plus = 0.5 * (D1 + D0); + main0[start_idx1] = 1.0 - 2.0 * D_plus * nu0; - D_minus = 0.5 * (D_row[nx - 2] + D_row[0]); - main0[nx - 1] = 1.0 - ((D_minus * nu0)); + D_minus = 0.5 * (Dend1 + Dend); + main0[end_idx1] = 1.0 - 2.0 * D_minus * nu0; - for (int j = 1; j < nx - 1; j++) + for (int j = start_idx2; j < end_idx1; j++) { D_minus = 0.5 * (D_row[j - 1] + D_row[j]); D_plus = 0.5 * (D_row[j + 1] + D_row[j]); main0[j] = 1.0 - ((D_minus * nu0) + (D_plus * nu0)); } - for (int j = 1; j < nx; j++) + for (int j = start_idx2; j < end_idx1+1; j++) { D_minus = 0.5 * ( D_row[j] + D_row[j - 1]); off_d_val_l0[j - 1] = nu0 * D_minus; } - for (int j = 0; j < nx - 1; j++) + for (int j = start_idx1; j < end_idx1; j++) { D_plus = 0.5 * (D_row[j + 1] + D_row[j]); off_d_val_u0[j] = nu0 * D_plus; } A[i] = new TridiagonalMatrix(main0, off_d_val_l0, off_d_val_u0); + #endregion ////================================================ //// The B-row Matrix next ////================================================ - //main = 1 + (2 * nu); - - D_plus = 0.5 * (D_row[1] + D_row[0]); - main1[0] = 1.0 + ((D_plus * nu0)); + #region B_row + D_plus = 0.5 * (D1 + D0); + main1[start_idx1] = 1.0 + 2.0 * D_plus * nu0; - D_minus = 0.5 * (D_row[nx - 2] + D_row[0]); - main1[nx - 1] = 1.0 + ((D_minus * nu0)); + D_minus = 0.5 * (Dend1 + Dend); + main1[end_idx1] = 1.0 + 2.0 * D_minus * nu0; - for (int j = 1; j < nx - 1; j++) + for (int j = start_idx2; j < end_idx1; j++) { D_minus = 0.5 * (D_row[j - 1] + D_row[j]); D_plus = 0.5 * (D_row[j + 1] + D_row[j]); main1[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0)); } - for (int j = 1; j < nx; j++) + for (int j = start_idx2; j < end_idx1+1; j++) { D_minus = 0.5 * (D_row[j - 1] + D_row[j]); off_d_val_l1[j - 1] = -nu0 * D_minus; } - for (int j = 0; j < nx - 1; j++) + for (int j = start_idx1; j < end_idx1; j++) { D_plus = 0.5 * (D_row[j + 1] + D_row[j]); off_d_val_u1[j] = -nu0 * D_plus; } B_row[i] = new TridiagonalMatrix(main1, off_d_val_l1, off_d_val_u1); + #endregion ////================================================ //// The B-column Matrix last ////================================================ + #region B_col RVector D_col = D.GetColVector(i); - - D_plus = 0.5 * (D_col[1] + D_col[0]); - main2[0] = 1.0 + ((D_plus * nu0)); + D0 = D_col[0]; + D1 = D_col[1]; + Dend = D_col[nx - 1]; + Dend1 = D_col[nx - 2]; - D_minus = 0.5 * (D_col[nx - 2] + D_col[0]); - main2[nx - 1] = 1.0 + ((D_minus * nu0)); + D_plus = 0.5 * (D1 + D0); + main2[start_idx1] = 1.0 + 2.0 * D_plus * nu0; - for (int j = 1; j < nx - 1; j++) + D_minus = 0.5 * (Dend1 + Dend); + main2[end_idx1] = 1.0 + 2.0 * D_minus * nu0; + + for (int j = start_idx2; j < end_idx1; j++) { D_minus = 0.5 * (D_col[j - 1] + D_col[j]); D_plus = 0.5 * (D_col[j + 1] + D_col[j]); main2[j] = 1.0 + ((D_minus * nu0) + (D_plus * nu0)); } - for (int j = 1; j < nx; j++) + for (int j = start_idx2; j < end_idx1+1; j++) { D_minus = 0.5 * (D_col[j - 1] + D_col[j]); off_d_val_l2[j - 1] = -nu0 * D_minus; } - for (int j = 0; j < nx - 1; j++) + for (int j = start_idx1; j < end_idx1; j++) { D_plus = 0.5 * (D_col[j + 1] + D_col[j]); off_d_val_u2[j] = -nu0 * D_plus; } - ////for (int j = 0; j < nx - 1; j++) { off_d_val_l[j] = nu[j]; } - ////for (int j = 1; j < nx; j++) { off_d_val_u[j - 1] = nu[j]; } - B_col[i] = new TridiagonalMatrix(main2, off_d_val_l2, off_d_val_u2); + #endregion } + // ==================================== // Set-up the initial condition function and start establishing the initial composition field // =================================== @@ -238,10 +258,10 @@ public DiffusionSimulators_2D_MatrixD( { border_with_function[i].PositionVaries = X.GetRowVector(X.GetnRows - 1); border_with_function[i].PositionFixed = Y[X.GetnRows - 1, 0]; - //C0 = border_with_function[i].BoundaryFunction(0.0, border_with_function[i].PositionVaries, border_with_function[i].PositionFixed); + C0 = border_with_function[i].BoundaryFunction(0.0, border_with_function[i].PositionVaries, border_with_function[i].PositionFixed); - //RVector Ctab = C_Initial.GetRowVector(ny - 1) + C0; - //C_Initial.ReplaceRow(C0, ny - 1); + RVector Ctab = C_Initial.GetRowVector(ny - 1) + C0; + C_Initial.ReplaceRow(C0, ny - 1); //Console.WriteLine(C0[0].ToString() + "Neumann Top!"); } break; @@ -402,8 +422,8 @@ public void Solve(int n_time_steps, int output_interval) cUpper2 = C_Im2[nrows - 2, ncols - 2]; } CR = RVector.Product(-2 * nu0 * cf_2D.DValues.GetColVector(ncols - 1), C1) + RVector.Product(2 * nu0 * cf_2D.DValues.GetColVector(ncols - 2), C2); // + C3 - CR[0] = ((-2 * nu0 * cf_2D.DValues[0, ncols - 1]) * cLower1) + ((2 * nu0 * cf_2D.DValues[1, ncols - 2]) * cLower2); - CR[ncols - 1] = ((-2 * nu0 * cf_2D.DValues[nrows - 1, ncols - 1]) * cUpper1) + ((2 * nu0 * cf_2D.DValues[nrows - 2, ncols - 2]) * cUpper2); + CR[0] = (-2 * nu0 * cf_2D.DValues[0, ncols - 1] * cLower1) + (2 * nu0 * cf_2D.DValues[1, ncols - 2] * cLower2); + CR[ncols - 1] = (-2 * nu0 * cf_2D.DValues[nrows - 1, ncols - 1] * cUpper1) + (2 * nu0 * cf_2D.DValues[nrows - 2, ncols - 2] * cUpper2); break; default: break; @@ -437,8 +457,8 @@ public void Solve(int n_time_steps, int output_interval) cUpper2 = C_Im2[nrows - 2, 1]; } CL = RVector.Product(-2 * nu0 * cf_2D.DValues.GetColVector(0), C1) + RVector.Product(2 * nu0 * cf_2D.DValues.GetColVector(1), C2); // + C3 - CL[0] = ((-2 * nu0 * cf_2D.DValues[0, 0]) * cLower1) + ((2 * nu0 * cf_2D.DValues[1, 1]) * cLower2); - CL[ncols - 1] = ((-2 * nu0 * cf_2D.DValues[nrows - 1, 0]) * cUpper1) + ((2 * nu0 * cf_2D.DValues[nrows - 2, 1]) * cUpper2); + CL[0] = (-2 * nu0 * cf_2D.DValues[0, 0] * cLower1) + (2 * nu0 * cf_2D.DValues[1, 1] * cLower2); + CL[ncols - 1] = (-2 * nu0 * cf_2D.DValues[nrows - 1, 0] * cUpper1) + (2 * nu0 * cf_2D.DValues[nrows - 2, 1] * cUpper2); break; default: @@ -475,9 +495,6 @@ public void Solve(int n_time_steps, int output_interval) C_Ex[i, 0] = CL[i]; //nu * C_Ex[i, ncols - 1] = CR[i]; //nu * } - //Console.WriteLine(t.ToString() + "s"); - //Console.WriteLine(BCs_Functions[0].BoundaryFunction(t * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed)[0].ToString()); - //Console.ReadKey(); // =================== // =================== // One-half implicit time-step @@ -487,7 +504,7 @@ public void Solve(int n_time_steps, int output_interval) if (t == 0) { fn = gxt_function(X, Y, t12, cf_2D.DValues, C_Initial); } else { fn = gxt_function(X, Y, t12, cf_2D.DValues, C_Im2); } - f12 = (dt / 2.0) * fn; + f12 = dt / 2.0 * fn; // BCs switch (BCs_Functions[0].TypeBC) @@ -495,7 +512,7 @@ public void Solve(int n_time_steps, int output_interval) case ABoundaryCondition.dirichlet: RVector gn = BCs_Functions[0].BoundaryFunction((t + 1) * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed); RVector g0 = BCs_Functions[0].BoundaryFunction(t * dt, BCs_Functions[0].PositionVaries, BCs_Functions[0].PositionFixed); - CT = ((0.5 * B_row[nrows - 1].Dot(gn)) + (0.5 * A[nrows - 1].Dot(g0))); + CT = (0.5 * B_row[nrows - 1].Dot(gn)) + (0.5 * A[nrows - 1].Dot(g0)); break; case ABoundaryCondition.neumann: RVector C1, C2; //,C3; @@ -532,7 +549,7 @@ public void Solve(int n_time_steps, int output_interval) case ABoundaryCondition.dirichlet: RVector gn = BCs_Functions[3].BoundaryFunction((t + 1) * dt, BCs_Functions[3].PositionVaries, BCs_Functions[3].PositionFixed); RVector g0 = BCs_Functions[3].BoundaryFunction(t * dt, BCs_Functions[3].PositionVaries, BCs_Functions[3].PositionFixed); - CB = ((0.5 * B_row[0].Dot(gn)) + (0.5 * A[0].Dot(g0))); //((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu * + CB = (0.5 * B_row[0].Dot(gn)) + (0.5 * A[0].Dot(g0)); //((0.5 * Bm.Dot(gn)) + (0.5 * Bp.Dot(g0))); //nu * break; case ABoundaryCondition.neumann: RVector C1, C2; //,C3; diff --git a/RMatrix.cs b/RMatrix.cs index 5c4cf3b..79a6621 100644 --- a/RMatrix.cs +++ b/RMatrix.cs @@ -642,7 +642,7 @@ public static RMatrix Adjoint(RMatrix rm) { for (int j = 0; j < rm.GetnCols; j++) { - ma[i, j] = Math.Pow(-1, i + j) * (Determinant(Minor(rm, i, j))); + ma[i, j] = Math.Pow(-1, i + j) * Determinant(Minor(rm, i, j)); } } return ma.GetTranspose(); @@ -656,7 +656,7 @@ public static RMatrix Inverse(RMatrix rm) { throw new Exception("Cannot inverse a matrix with 0 determinant!"); } - return (Adjoint(rm) / Dm); + return Adjoint(rm) / Dm; } public static RMatrix InverseAccurate(RMatrix rm) { @@ -811,7 +811,7 @@ public static void Cholesky_Decomposition(RMatrix rm) { // Evaluating L(i, j) // using L(j, j) - for (int k = 0; k < j; k++) { sum += (lower[i, k] * lower[j, k]); } + for (int k = 0; k < j; k++) { sum += lower[i, k] * lower[j, k]; } lower[i, j] = (rm[i, j] - sum) / lower[j, j]; } } @@ -910,7 +910,7 @@ public static RMatrix GaussSeidel(RMatrix A, RVector x0, RMatrix b) sum = b[i, 0]; for (int j = 0; j < ncols; j++) { - if (j != i) { sum -= (A[i, j] * x[j]); } + if (j != i) { sum -= A[i, j] * x[j]; } } x[i] = sum / A[i, i]; } diff --git a/TridiagonalMatrix.cs b/TridiagonalMatrix.cs index eaa059a..922c081 100644 --- a/TridiagonalMatrix.cs +++ b/TridiagonalMatrix.cs @@ -222,7 +222,7 @@ public static RVector GaussSeidel_Old(TridiagonalMatrix A, RVector xold, RVector { Rval += A[i, k] * xold[k]; } - xnew[i] = (1.0 / A[i, i]) * (b[i] - Lval - Rval); + xnew[i] = 1.0 / A[i, i] * (b[i] - Lval - Rval); } return xnew; @@ -280,7 +280,7 @@ public static RVector GaussSeidel(TridiagonalMatrix A, RVector b) //RVector XO, sum = b[j] - (Ajjm1 * x[j - 1]) - (Ajjp1 * XO[j + 1]); } Ajj = A[j, j]; - x[j] = (1.0 / Ajj) * (sum); + x[j] = 1.0 / Ajj * sum; } RVector normx = x - XO; double normx_val = normx.GetNorm(); @@ -319,7 +319,7 @@ public static RVector SOR(TridiagonalMatrix A, RVector b) sum = b[j] - (Ajjm1 * x[j - 1]) - (Ajjp1 * XO[j + 1]); } Ajj = A[j, j]; - x[j] = ((1 - omega) * XO[j]) + (omega / Ajj) * (sum); + x[j] = ((1 - omega) * XO[j]) + omega / Ajj * sum; } RVector normx = x - XO; double normx_val = normx.GetNorm(); diff --git a/log.txt b/log.txt index f7254a6..9089518 100644 --- a/log.txt +++ b/log.txt @@ -242,3 +242,17 @@ {"message":"No project detected for extracting metadata.","source":"MetadataCommand.ExtractMetadata","date_time":"2021-07-21T19:01:00.9069456Z","message_severity":"warning","correlation_id":"DDA2B50C-E1AD-49DC-A6EB-33DDC392F651.1.1.3"} {"message":"Unable to find either toc.yml or toc.md inside api/. Make sure the file is included in config file docfx.json!","source":"BuildCommand.BuildCore.Build Document.CompilePhaseHandlerWithIncremental.TocDocumentProcessor.Prebuild.BuildTocDocument","file":"toc.yml","date_time":"2021-07-21T19:01:01.90895Z","message_severity":"warning","correlation_id":"DDA2B50C-E1AD-49DC-A6EB-33DDC392F651.2.180.1.31.8.3.2.1"} {"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-21T19:01:02.3309562Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T15:03:25.3848108Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T15:32:07.1473506Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T16:36:42.0548635Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T16:37:36.0979328Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T16:40:59.4669517Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T16:44:03.4841673Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T16:44:54.6011534Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T16:56:41.6962107Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T17:15:59.909172Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T17:22:37.2516399Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T17:24:29.3370322Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T17:26:21.8812888Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T17:47:43.7423873Z","message_severity":"warning","code":"InvalidFileLink"} +{"message":"Invalid file link:(~/2D_Region.png).","source":"BuildCommand.BuildCore.Build Document.LinkPhaseHandlerWithIncremental.ConceptualDocumentProcessor.Save","file":"README.md","line":"0","date_time":"2021-07-22T17:53:02.6368034Z","message_severity":"warning","code":"InvalidFileLink"}