Clear["Global`*"];
delta[a_] := Z^2*(1 - 2 a)/(a*(Z + 1));
a11[a_] := (3*Pi/32)*(288 + 604*Sqrt[2]*delta[a] + 217*delta[a]^2)/
4/(72 + 61*Sqrt[2]*delta[a] + 16*delta[a]^2);
dc[a_] :=
Piecewise[{{1,
a <= 0}, {a11[a]/(1 - a*(2 - mi/mI*(Z + 1)))/(a*(Z - 1) + 1),
a > 0 && a <= 1/2}, {dc[1/2], a > 1/2}}];
dp[a_] := dc[a]*a*(1 - 2 a)*(Z - 1)^2/(Z - a*(Z - 1));
D0 = 278*4/100*Sqrt[2]*5^(5/2)/(2*1.81)*(Z + 1)/Z^2;
(*cm^2/s*)
Z = 6;
mi = 2; mI = 36;
a0[z_] := Piecewise[{{1/2, z < 0}, {0, z >= 0}}];
sol = NDSolve[{D[a[z, t], t] ==
D[D0*(dc[a[z, t]] + dp[a[z, t]])*D[a[z, t], z], z],
a[-100, t] == 1/2, a[100, t] == 0, a[z, 0] == a0[z]},
a, {z, -100, 100}, {t, 0, 3}]
这是一个扩散方程的数值求解。
在NDSolve的过程中，会遇到Divide:碰到无穷表达式1/0. 这个问题很显然来自dc[a]，当我将其设置为常数时可以得到很好的解。
具体原因应该是dc[a]中的delta[a]在a=0时分母为0。然而事实上，在[0,1/2]的定义域上，dc[a]的原始定义（也就是我分段函数中的中间那段）的表现都是良好的，画图也可以看出dc[0]=1;
我尝试将dc[a]设置为分段函数来解决这个问题，然而并没有作用。有什么好的办法来解决这个问题？我能想到把dc[a]手动插值一遍，但是既然dc[a]本身表现良好，我想应该有更好的办法。
delta[a_] := Z^2*(1 - 2 a)/(a*(Z + 1));
a11[a_] := (3*Pi/32)*(288 + 604*Sqrt[2]*delta[a] + 217*delta[a]^2)/
4/(72 + 61*Sqrt[2]*delta[a] + 16*delta[a]^2);
dc[a_] :=
Piecewise[{{1,
a <= 0}, {a11[a]/(1 - a*(2 - mi/mI*(Z + 1)))/(a*(Z - 1) + 1),
a > 0 && a <= 1/2}, {dc[1/2], a > 1/2}}];
dp[a_] := dc[a]*a*(1 - 2 a)*(Z - 1)^2/(Z - a*(Z - 1));
D0 = 278*4/100*Sqrt[2]*5^(5/2)/(2*1.81)*(Z + 1)/Z^2;
(*cm^2/s*)
Z = 6;
mi = 2; mI = 36;
a0[z_] := Piecewise[{{1/2, z < 0}, {0, z >= 0}}];
sol = NDSolve[{D[a[z, t], t] ==
D[D0*(dc[a[z, t]] + dp[a[z, t]])*D[a[z, t], z], z],
a[-100, t] == 1/2, a[100, t] == 0, a[z, 0] == a0[z]},
a, {z, -100, 100}, {t, 0, 3}]
这是一个扩散方程的数值求解。
在NDSolve的过程中，会遇到Divide:碰到无穷表达式1/0. 这个问题很显然来自dc[a]，当我将其设置为常数时可以得到很好的解。
具体原因应该是dc[a]中的delta[a]在a=0时分母为0。然而事实上，在[0,1/2]的定义域上，dc[a]的原始定义（也就是我分段函数中的中间那段）的表现都是良好的，画图也可以看出dc[0]=1;
我尝试将dc[a]设置为分段函数来解决这个问题，然而并没有作用。有什么好的办法来解决这个问题？我能想到把dc[a]手动插值一遍，但是既然dc[a]本身表现良好，我想应该有更好的办法。
