I am trying to numerically solve an integral equation in Starfield and Crouch textbook (Boundary Element Methods in Solid Mechanics: With Applications in Rock Mechanics and Geological Engineering), equation 3.2.5 (see code). After computing the density py
I attempted to plot uy
vs x
. I observed some numerical instabilities in the region 0<=x<=1
. I Increased my number of integration points,np
and also reduced my boundary element size,xe
, I still observe the same irregularity.I'm not sure why this is happening. Is there something I'm not doing right? I would appreciate any guidance. Here is my code:
(I have added the corresponding analytic solution plot)
(*Boundary elements set up and material properties*)
nb = nm = 20; nd=ne=nb+1; G = 1; v = 0.1; L = 1.25;a = 0.1;
(*Input the coordinates of the of ends of boundary elements (xe,ye)*)
xe = Table[i, {i, -1, 1, 2/nb}];
ye = Table[0, {i, 1, ne}];
(*Input the coordinates of the midpoints of boundary elements (xm,ym)*)
xm = ym = Table[0, {i, 1, nm}];
jb = If[j < ne, j + 1];
Do[xm[[j]] = (xe[[j]] + xe[[jb]])/2;
ym[[j]] = (ye[[j]] + ye[[jb]])/2, {j, 1, nm}];
bv = Table[-1, {i, 1, nm}];
(*Compute elements of Influence coefficients Bij and Sij*)
Sij = Bij = Table[0, {i, 1, nb}, {j, 1, nb}];
uy = (1/(2 G Pi)) (-2 (1 - v) (Log[Sqrt[(x - xi)^2 + y^2]] -
Log[L - xi]) + y^2/((x - xi)^2 + y^2))(*Equation 3.2.5*);
Get["NumericalDifferentialEquationAnalysis`"];
np = 6; points = weights = Table[Null, {np}];
Do[points[[i]] = GaussianQuadratureWeights[np, -1, 1][[i, 1]], {i, 1, np}]
Do[weights[[i]] = GaussianQuadratureWeights[np, -1, 1][[i, 2]], {i, 1, np}]
GuassInt[f_, z_] := Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}]
Do[xb = (1/2)*(xe[[jb]]*(1 - z) + xe[[j]]*(1 + z)); yb = (1/2)*(ye[[jb]]*(1 -z) + ye[[j]]*(1 + z));
Do[Bij[[i, j]] = GuassInt[uy /. {x -> xm[[i]], xi -> xb, y -> yb}, z]; Sij[[i, j]] = GuassInt[uy /. {x -> x, y -> yb, xi -> xb}, z], {i, 1, nb}], {j, 1, nb}]
py = LinearSolve[Bij, bv];
plot1 = Plot[Sij . py, {x, 0, 3},PlotStyle -> Blue]
AnalyticUy[h_] := -(1 - ((Log[h + Sqrt[(h^2) - 1]])/Log[2]))
plot2 = Plot[AnalyticUy[h], {h, 1, 3}, PlotStyle -> {Dashed, Black}]
plot3 = Plot[-1, {x, 0, 1}, PlotStyle -> {Dashed, Black}]
Show[plot1, plot2, plot3]
np
to20
, I find the irregularity has reduced significantly. $\endgroup$ – xzczd Dec 23 '17 at 4:40xe
has the wrong number of elements. Probably you would like to havexe = Table[i, {i, -1, 1, 2/nb}];
in your code. Moreover, I would advise you to computend
,ne
andnm
fromnb
if possible; this way you have fewer place to change whenever you change parameters. $\endgroup$ – Henrik Schumacher Dec 23 '17 at 23:25