NDSolve 无法求解偏微分方程

Question

求解偏微分方程时，我得到

NDSolve::eerri：警告：在自变量 x 方向上指定空间网格上的估计初始误差超出了规定的误差容限。并且方程无法解。

我的代码：

ClearAll
mu = 0.25
rho = 2350
xE = 18*1000000000
G = xE/(2*(1 + mu))
cp = Sqrt[xE*(1 - mu)/(rho*(1 + mu)*(1 - 2*mu))]
r0 = 2.5
xdb = 40/1000
xdp = 2.5*40/1000
rho0 = 1100
xD = 3200
sb = 4500/1000
r0 = 2.5
b = 1.5
alpha = 4000
JM = (3*Pi)/(4*r0)
m = alpha/cp
j = Sqrt[JM^2 + m^2]
A0 = xdb/(8*sb)*rho0*xD^2/(xdp/xdb)^2.2
A = (1 - E^(-b*j))*
  xE ((j^2 - m^2 + (2*j^2 - 3*m^2)*mu)*r0^(-1/2)*
      Cos[Sqrt[j^2 - m^2]*r0 - Pi/4]/(1 - 2*mu) - 
     Sqrt[j^2 - m^2]*r0^(-3/2)*
      Sin[Sqrt[j^2 - m^2]*r0 - Pi/4] - (3 - 5*mu)*r0^(-5/2)*
      Cos[Sqrt[j^2 - m^2]*r0 - Pi/4]/(4*(1 - 2*mu)))/(b*j*(1 + mu))
A1 = A0/A
Ur[r_, z_, t_] = 
 FullSimplify[
  A1*(1/2*Cos[Sqrt[j^2 - m^2]*r - 1/4*Pi]/r^(3/2) + 
     Sqrt[j^2 - m^2]*Sin[Sqrt[j^2 - m^2]*r - 1/4*Pi]/r^(1/2))*
   E^(-j*z)*E^(-m*cp*t)]
Vr[r_, z_, t_] = 
 FullSimplify[-m*cp*
   A1*(1/2*Cos[Sqrt[j^2 - m^2]*r - 1/4*Pi]/r^(3/2) + 
     Sqrt[j^2 - m^2]*Sin[Sqrt[j^2 - m^2]*r - 1/4*Pi]/r^(1/2))*
   E^(-j*z)*E^(-m*cp*t)]
Sigmaz[r_, z_, t_] = 
 FullSimplify[
  xE*A1*((-j^2 + (2*j^2 - m^2)*mu)*
      Cos[Sqrt[j^2 - m^2]*r - 1/4*Pi]/((1 - 2*mu)*r^(1/2)) - 
     mu*Cos[Sqrt[j^2 - m^2]*r - 1/4*Pi]/((4 - 8*mu)*r^(5/2)))*
   E^(-j*z)*E^(-m*cp*t)/(1 + mu)]
l = 2.5
l0 = 1.5
Eb = 210*1000000000
rb = 10/1000
rhob = 7850
Eg = 35*1000000000
rg = 17.5/1000
vg = 0.25
varphi = 45/180*Pi
Cb = 2*Pi*rb
Ab = Pi*rb^2
R = 10*Eb*rb/((Eg + xE)/2)
Gg = Eg/(2*(1 + vg))
k = Gg*G/((G*Log[rg/rb] + Gg*Log[R/rb])*rb)
YiTa = Tan[varphi]
z1 = b
p[x_, t_] = FullSimplify[YiTa*Sigmaz[r0 + l - x, z1, t]]
PDE1 = D[u1[x, t], {x, 2}] == 
  rhob*D[u1[x, t], {t, 2}]/Eb + k*Cb*u1[x, t]/(Eb*Ab) + 
   p[x, t]*Cb/(Eb*Ab)
initialCondition = {Derivative[0, 1][u1][x, 0] == 
   Vr[r0 + l - x, z1, 0], u1[x, 0] == Ur[r0 + l - x, z1, 0]}
boundaryConditions = {u1[0, t] == Ur[l + r0, z1, t], 
  u1[l0, t] == Ur[r0 + l - l0, z1, t]}
sol = NDSolve[{PDE1, initialCondition, boundaryConditions}, 
  u1, {x, 0, l0}, {t, 0, 0.001}, MaxStepSize -> 0.0000001]
Plot3D[Evaluate[u1[x, t] /. sol], {x, 0, l0}, {t, 0, 0.0001}, 
 PlotRange -> All, AxesLabel -> {"x", "t", "u1(x,t)"}]

我尝试运行代码但没有结果。我希望得到可以运行和绘制的结果。

Answer 1

评论太长了。正如我之前所建议的，我用精确的有理数替换了您的

$MachinePrecision

常数。

mu=25/100;rho=2350;xE=18*10^9;G=xE/(2*(1+mu));cp=Sqrt[xE*(1-mu)/(rho*(1+mu)*(1-2*mu))];
r0=25/10;xdb=40/1000;xdp=25/10*40/1000;rho0=1100;xD=3200;sb=4500/1000;r0=25/10;
b=15/10;alpha=4000;
JM=(3*Pi)/(4*r0);m=alpha/cp;j=Sqrt[JM^2+m^2];A0=xdb/(8*sb)*rho0*xD^2/(xdp/xdb)^(22/10);
A=(1-E^(-b*j))*xE((j^2-m^2+(2*j^2-3*m^2)*mu)/Sqrt[r0]*Cos[Sqrt[j^2-m^2]*r0-
  Pi/4]/(1-2*mu)-Sqrt[j^2-m^2]*r0^(-3/2)*Sin[Sqrt[j^2-m^2]*r0-Pi/4]-(3-5*mu)*
  r0^(-5/2)*Cos[Sqrt[j^2-m^2]*r0-Pi/4]/(4*(1-2*mu)))/(b*j*(1+mu));A1=A0/A;
Ur[r_,z_,t_]=A1*(1/2*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/r^(3/2)+Sqrt[j^2-m^2]*
  Sin[Sqrt[j^2-m^2]*r-1/4*Pi]/Sqrt[r])*E^(-j*z)*E^(-m*cp*t);
Vr[r_,z_,t_]:=-m*cp*A1*(1/2*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/r^(3/2)+Sqrt[j^2-
  m^2]*Sin[Sqrt[j^2-m^2]*r-1/4*Pi]/Sqrt[r])*E^(-j*z)*E^(-m*cp*t);
Sigmaz[r_,z_,t_]:=xE*A1*((-j^2+(2*j^2-m^2)*mu)*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/
  ((1-2*mu)*Sqrt[r])-mu*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/((4-8*mu)*r^(5/2)))*
  E^(-j*z)*E^(-m*cp*t)/(1+mu);
l=25/10;l0=15/10;Eb=210*10^9;rb=10/1000;rhob=7850;Eg=35*10^9;rg=175/10/1000;
vg=025/100;varphi=45/180*Pi;Cb=2*Pi*rb;Ab=Pi*rb^2;R=10*Eb*rb/((Eg+xE)/2);
Gg=Eg/(2*(1+vg));k=Gg*G/((G*Log[rg/rb]+Gg*Log[R/rb])*rb);YiTa=Tan[varphi];z1=b;
p[x_,t_]:=YiTa*Sigmaz[r0+l-x,z1,t];
PDE1={D[u1[x,t],{x,2}]==rhob*D[u1[x,t],t,2}]/Eb+k*Cb*u1[x,t]/(Eb*Ab)+
  p[x,t]*Cb/(Eb*Ab)};
initialCondition={Derivative[0,1][u1][x,0]==Vr[r0+l-x,z1,0],u1[x,0]==Ur[r0+l-x,z1,0]};
boundaryConditions={u1[0,t]==Ur[l+r0,z1,t],u1[l0,t]==Ur[r0+l-l0,z1,t]};
sys=Simplify[N[Join[PDE1,initialCondition,boundaryConditions],64]]

就在执行（大致）显示的

NDSolve

的步骤之前

{Derivative[2, 0][u1][x, t] == (     
0.00545860389386734035988805739383977751117746615039472979478311832031082160733*
(25.035720659226906140411453104565220866207337272409961236398875988369353985859-
10.*x + 1.*x^2)*Cos[
3.9269908169872415480783042290993786052464617492188822762186807403847705078577- 
.94247779607693797153879301498385086525915081981253174629248337769234492*x])/    
(E^(3999.99999999999999999999999999999999999999999999999999999999999999999*t)* 
(5. - 1.*x)^(5/2)) +    
147.137927874474292837143887032039529266866273227161600896429319953505602294*
u1[x, t] + 3.738095238095238095238095238095238095238095238095238095238095238*^
-8*Derivative[0, 2][u1][x, t], 
Derivative[0, 1][u1][x, 0] == -(((       
(1.6768709894289540203312443887507074674202073320420657991243447623718564+ 0.*I)- 
0.30320332583960679232662864081734027262650464178210400129298177083547+0.*I)*x)*
Cos[0.9424777960769379715387930149838508652591508198125317462924833776923*x] + 
((-1.3551622689671139029350420194226952588448390857789742138054729459828 +0.*I) + 
(0.303203325839606792326628640817340272626504641782104001292981770835+0.*I)*x)*
Sin[0.94247779607693797153879301498385086525915081981253174629248337769*x])/
(5. - 1.*x)^(3/2)), 
u1[x, 0] == -((( 
(-0.00041921774735723850508281109718767686685505183301051644978108619059+0.*I) +
(0.00007580083145990169808165716020433506815662616044552600032324544 +0.*I)*x)*
Cos[0.9424777960769379715387930149838508652591508198125317462924833776923*x] + 
((0.000338790567241778475733760504855673814711209771444743553451368236495+0.*I)- 
(0.0000758008314599016980816571602043350681566261604455260003232454427+0.*I)*x)*
Sin[0.942477796076937971538793014983850865259150819812531746292483377688587*x])/
(5. - 1.*x)^(3/2)), 
u1[0, t] == 
0.000037495975218604724441921501083994167027199906049745772524411502557/
E^(3999.99999999999999999999999999999999999999999999999999999999999999999993*t),
u1[1.5, t] == 
-0.00002665354928541077808491725593466020104224160610364654413085048230532/
E^(3999.9999999999999999999999999999999999999999999999999999999999999999993*t)}

除了 E^(3999.99999...*t) 项外，大部分看起来都是合理的 以及所有“大约 0*I”项。

你能完成我所做的更改，看看我这样做是否犯了任何错误，并删除那些让我害怕的术语吗？

我还没有删除这些，但我尝试了过去使用过的一些技巧，它要么花费大量时间或内存，要么通过有关错误太大的警告来解决。