背景:我正在尝试求解一个不稳定性问题,它最终被我归结为含有一个未知参数w(也是我希望求解的值)的两个耦合微分方程组。由于这个问题的复杂性,许多传统方法(Shooting method,本征矩阵法)都不太合适。我目前想到的最好方法是这样的:
ClearAll["Global`*"];
k = Sqrt[3/8]*wp/V0; wp = 1; wc = wp*V0/c;
c = 1; V0 = 1/5; gamma0 := 1/Sqrt[1 - V0^2/c^2];
w0[k_] :=
N[Sqrt[wp^2/2/
gamma0^3*(Sqrt[1 + 8*k^2*V0^2/wp^2*gamma0^3] - 1 -
2*k^2*V0^2/wp^2*gamma0^3)]];
Omega[V_] := w - k*V;
A[V_] := (w - k*V)^2*(w^2 - k^2*c^2 - wp^2/gamma0) -
wc^2/gamma0^2*(w^2 - k^2*c^2);
B[V_] :=
1/c^2*((w - k*V)^2*(w^2 - k^2*c^2 - wp^2/gamma0)^2 -
wc^2/gamma0^2*(w^2 - k^2*c^2)^2);
Eq1u = w^2/A[V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[V0]^2)*
Eyu''[x] +
w^2/c^2 (wp^2/gamma0^3/Omega[V0]^2 - 1)*
Eyu[x] - (w*wc/gamma0*wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*
Ezu'[x] == 0;
Eq1d = w^2/A[-V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[-V0]^2)*
Eyd''[x] +
w^2/c^2 (wp^2/gamma0^3/Omega[-V0]^2 - 1)*
Eyd[x] - (w*wc/gamma0*wp^2/gamma0)*(k + w*V0/c^2)/A[-V0]*
Ezd'[x] == 0;
Eq2u = Ezu''[x] +
B[V0]/A[V0]*
Ezu[x] - (w*wc/gamma0*wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*
Eyu'[x] == 0;
Eq2d = Ezd''[x] +
B[V0]/A[V0]*
Ezd[x] - (w*wc/gamma0*wp^2/gamma0)*(k + w*V0/c^2)/A[V0]*
Eyd'[x] == 0;
Eq1c = w^2/A[V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[V0]^2)*
Eyu'[0] - (w*wc/gamma0*wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*Ezu[0] ==
w^2/A[-V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[-V0]^2)*
Eyd'[0] - (w*wc/gamma0*wp^2/gamma0)*(k + w*V0/c^2)/A[-V0]*
Ezd[0];
Eq2c = Ezu'[0] == Ezd'[0];
yc = Eyu[0] == Eyd[0];
zc = Ezu[0] == Ezd[0];
sol = N[DSolve[{Eq1u, Eq2u, Eq1d, Eq2d}, {Eyu, Ezu, Eyd, Ezd}, x]];
eqns = {First[Ezu[0] /. sol] == First[Ezd[0] /. sol],
First[Ezu'[0] /. sol] == First[Ezd'[0] /. sol],
First[Eyu[0] /. sol] == First[Eyd[0] /. sol],
w^2/A[V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[V0]^2)*
First[Eyu'[0] /. sol] - (w*wc/gamma0*
wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*First[Ezu[0] /. sol] ==
w^2/A[-V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[-V0]^2)*
First[Eyd'[0] /. sol] - (w*wc/gamma0*
wp^2/gamma0)*(k + w*V0/c^2)/A[-V0]*First[Ezd[0] /. sol],
C[4] == 1,
First[Eyu[Infinity] /. sol] == 0, First[Ezu[Infinity] /. sol] == 0,
First[Eyd[-Infinity] /. sol] == 0, First[Ezd[-Infinity] /. sol] == 0};
NSolve[eqns]
其中最后几个方程来自于不稳定模式在无穷远处消失的条件。但问题似乎也出在这四个方程上,因为w是未知的。我也尝试了把这四个方程中的Infinity换成一个较大的数值(例如100),但是这似乎会因为较大的系数造成问题。
在一些传统的问题里,不稳定模式在无穷远处消失的条件被用来消除形如c*Exp[a*x](a>0)的解。事实上,如果我代入一个w的试探值,例如
Simplify[Eyu[x] /. sol /. w -> I*w0[k]]
Out[97]= {E^((3.15224 + 0.0994554 I) x) ((0.209468 + 0.121609 I) C[
1] + (0.0651416 + 0.0441776 I) C[2] - (1.54602 - 1.42907 I) C[
3] - (0.474889 - 0.417018 I) C[4]) +
E^((-3.15224 - 0.0994554 I) x) ((0.209468 + 0.121609 I) C[
1] - (0.0651416 + 0.0441776 I) C[2] + (1.54602 - 1.42907 I) C[
3] - (0.474889 - 0.417018 I) C[4]) +
E^((-3.19809 + 0.152718 I) x) ((0.290532 - 0.121609 I) C[
1] - (0.095468 - 0.0410111 I) C[2] - (1.62954 - 1.28268 I) C[
3] + (0.474889 - 0.417018 I) C[4]) +
E^((3.19809 - 0.152718 I) x) ((0.290532 - 0.121609 I) C[
1] + (0.095468 - 0.0410111 I) C[2] + (1.62954 - 1.28268 I) C[
3] + (0.474889 - 0.417018 I) C[4])}
此时很明显,第一行和第四行的系数必须为0,这样就可以得到两个很简单的方程。
但是在w未知的情况下,我该怎么做这件事呢?
ClearAll["Global`*"];
k = Sqrt[3/8]*wp/V0; wp = 1; wc = wp*V0/c;
c = 1; V0 = 1/5; gamma0 := 1/Sqrt[1 - V0^2/c^2];
w0[k_] :=
N[Sqrt[wp^2/2/
gamma0^3*(Sqrt[1 + 8*k^2*V0^2/wp^2*gamma0^3] - 1 -
2*k^2*V0^2/wp^2*gamma0^3)]];
Omega[V_] := w - k*V;
A[V_] := (w - k*V)^2*(w^2 - k^2*c^2 - wp^2/gamma0) -
wc^2/gamma0^2*(w^2 - k^2*c^2);
B[V_] :=
1/c^2*((w - k*V)^2*(w^2 - k^2*c^2 - wp^2/gamma0)^2 -
wc^2/gamma0^2*(w^2 - k^2*c^2)^2);
Eq1u = w^2/A[V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[V0]^2)*
Eyu''[x] +
w^2/c^2 (wp^2/gamma0^3/Omega[V0]^2 - 1)*
Eyu[x] - (w*wc/gamma0*wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*
Ezu'[x] == 0;
Eq1d = w^2/A[-V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[-V0]^2)*
Eyd''[x] +
w^2/c^2 (wp^2/gamma0^3/Omega[-V0]^2 - 1)*
Eyd[x] - (w*wc/gamma0*wp^2/gamma0)*(k + w*V0/c^2)/A[-V0]*
Ezd'[x] == 0;
Eq2u = Ezu''[x] +
B[V0]/A[V0]*
Ezu[x] - (w*wc/gamma0*wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*
Eyu'[x] == 0;
Eq2d = Ezd''[x] +
B[V0]/A[V0]*
Ezd[x] - (w*wc/gamma0*wp^2/gamma0)*(k + w*V0/c^2)/A[V0]*
Eyd'[x] == 0;
Eq1c = w^2/A[V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[V0]^2)*
Eyu'[0] - (w*wc/gamma0*wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*Ezu[0] ==
w^2/A[-V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[-V0]^2)*
Eyd'[0] - (w*wc/gamma0*wp^2/gamma0)*(k + w*V0/c^2)/A[-V0]*
Ezd[0];
Eq2c = Ezu'[0] == Ezd'[0];
yc = Eyu[0] == Eyd[0];
zc = Ezu[0] == Ezd[0];
sol = N[DSolve[{Eq1u, Eq2u, Eq1d, Eq2d}, {Eyu, Ezu, Eyd, Ezd}, x]];
eqns = {First[Ezu[0] /. sol] == First[Ezd[0] /. sol],
First[Ezu'[0] /. sol] == First[Ezd'[0] /. sol],
First[Eyu[0] /. sol] == First[Eyd[0] /. sol],
w^2/A[V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[V0]^2)*
First[Eyu'[0] /. sol] - (w*wc/gamma0*
wp^2/gamma0)*(k - w*V0/c^2)/A[V0]*First[Ezu[0] /. sol] ==
w^2/A[-V0]*(wc^2/gamma0^2 + wp^2/gamma0^3 - Omega[-V0]^2)*
First[Eyd'[0] /. sol] - (w*wc/gamma0*
wp^2/gamma0)*(k + w*V0/c^2)/A[-V0]*First[Ezd[0] /. sol],
C[4] == 1,
First[Eyu[Infinity] /. sol] == 0, First[Ezu[Infinity] /. sol] == 0,
First[Eyd[-Infinity] /. sol] == 0, First[Ezd[-Infinity] /. sol] == 0};
NSolve[eqns]
其中最后几个方程来自于不稳定模式在无穷远处消失的条件。但问题似乎也出在这四个方程上,因为w是未知的。我也尝试了把这四个方程中的Infinity换成一个较大的数值(例如100),但是这似乎会因为较大的系数造成问题。
在一些传统的问题里,不稳定模式在无穷远处消失的条件被用来消除形如c*Exp[a*x](a>0)的解。事实上,如果我代入一个w的试探值,例如
Simplify[Eyu[x] /. sol /. w -> I*w0[k]]
Out[97]= {E^((3.15224 + 0.0994554 I) x) ((0.209468 + 0.121609 I) C[
1] + (0.0651416 + 0.0441776 I) C[2] - (1.54602 - 1.42907 I) C[
3] - (0.474889 - 0.417018 I) C[4]) +
E^((-3.15224 - 0.0994554 I) x) ((0.209468 + 0.121609 I) C[
1] - (0.0651416 + 0.0441776 I) C[2] + (1.54602 - 1.42907 I) C[
3] - (0.474889 - 0.417018 I) C[4]) +
E^((-3.19809 + 0.152718 I) x) ((0.290532 - 0.121609 I) C[
1] - (0.095468 - 0.0410111 I) C[2] - (1.62954 - 1.28268 I) C[
3] + (0.474889 - 0.417018 I) C[4]) +
E^((3.19809 - 0.152718 I) x) ((0.290532 - 0.121609 I) C[
1] + (0.095468 - 0.0410111 I) C[2] + (1.62954 - 1.28268 I) C[
3] + (0.474889 - 0.417018 I) C[4])}
此时很明显,第一行和第四行的系数必须为0,这样就可以得到两个很简单的方程。
但是在w未知的情况下,我该怎么做这件事呢?