画一个三维的箭头，z方向的范围比x，y方向大很多，默认比例画出来的箭头形状是没问题的 Graphics3D[{Red, Arrowheads[0.1], Arrow[Tube[{{1, 1, -1}, {2, 2, 0}, {3, 3, -1}, {4, 4, 20}}, 0.05]]}] 但是不能展示比例这么奇怪的图，要把z方向压缩一下，这样就改变了箭头的形状，变得很奇怪 Graphics3D[{Red, Arrowheads[0.1], Arrow[Tube[{{1, 1, -1}, {2, 2, 0}, {3, 3, -1}, {4, 4, 20}}, 0.05]]}, BoxRatios -> {1, 1.5, 1}] 好像没有三维箭头设置长宽比的选项，只能改大小，请问怎么让下面的图片中箭头形状变

lz没有系统性的学过积分方程，想用mathematica求解积分方程，查阅帮助文档后发现，mathematica好像没有关于积分方程的文档，百度了一下，一维积分方程可以用Sovle求解，但对于我的问题却失效了，以下是代码： eqn = {Integrate[f[x]/(2\[CenterDot]\[Pi]) 1/(x1 - x), {x, 0, c}] == Subscript[V, \[Infinity]] \[Alpha], f[0] == Infinity, f[c] == 0} Solve[eqn, f[x]] 求问mathematica怎么求解积分方程qvq

代码是 Clear["Global`*"] m := 1.320; \[Sigma] := 0.200; \[Alpha]s := 0.35; \[Mu]R := m/4; \[Delta] := 10^-10; V[r_] := -4/3*\[Alpha]s/r + \[Sigma]*r; Veff[r_] := 2*\[Mu]R*V[r]; SolveRadialEq[\[CurlyEpsilon]_] := NDSolve[{u''[r] == (Veff[r] - \[CurlyEpsilon]^2/4 + m^2) *u[r], u[0] == \[Delta], u'[0] == 1}, u, {r, 0, 20}]; CheckBoundaryCondition[sol_, \[CurlyEpsilon]_] := Module[{uSol, value}, uSol := u[r] /. First[sol]; value := uSol[20]; value]; FindEnergyLevel[l_, \[CurlyEpsilon]Min_, \[CurlyEpsilon]Max_, tolerance_] := Module[{\[CurlyEpsilon]Mid, sol, check}, While[

此帖是本吧初代吧主 @mm_酱 发在果壳小组的一篇教程。如今果壳小组已死，mm_酱近期又不太活跃，这里姑且重发一份： https://note.youdao.com/share/?id=058e6037396d925af1f4abe4d54a52 mm_酱你要是想整理下出个v2就吱一声。

把求根公式换成chareqn2/chareqn3 跑的时候会显示计算溢出 搞不明白 有没有老哥救一下 底下放源码 Clear; \[Lambda] = 0.633; nSiO2 = Sqrt[ 1 + x1*\[Lambda]^2/(\[Lambda]^2 - x2^2) + x3*\[Lambda]^2/(\[Lambda]^2 - x4^2) + x5*\[Lambda]^2/(\[Lambda]^2 - x6^2) /. {x1 -> 0.6961663, x2 -> 0.0684043, x3 -> 0.4079426, x4 -> 0.1162414, x5 -> 0.8974794, x6 -> 9.896161}]; nSi = Sqrt[ 11.6858 + 0.939816/\[Lambda]^2 + 0.000993358/(\[Lambda]^2 - 1.22567)]; n1 = nSiO2; n2 = 1; k0 = 2*\[Pi]/\[Lambda] c = 2.99979*10^14; d0 = 0.12; (*0.12um=120nm*) \[CapitalDelta]d =

现在的情况是这样的，中间的错误不知道如何修改，请各位大佬能帮我看一下怎么修改，现在的情况是没办法出图，应该是中间的哪一步的参数代换错了？还是有程序哪一步没走通。 程序里应该是建立映射来管理特征值或者解的排序，不清楚，我这段映射的程序用的对不对。 麻烦各位大佬了。 p = K0^2 \[Omega] - Ke^2 \[Omega] + Kf^2 \[Omega] - Kv^2 \[Omega] + \[Epsilon] \[Omega]^2 - \[Omega]^3 - 2 Ke Kf Subscript[l, B] + 2 K0 Kv Subscript[l, B] + \[Epsilon] \!\(\*SubsuperscriptBox[\(l\), \(B\), \(2\

理论上来说,第一类曲面积分对1积分的结果应该是曲面的面积 然而今天计算x^2+y^2==Cos[z]^2的函数时,其积分结果总为0 同时尝试了一下x^2+y^2==z^2,积分结果是正常的 然后又尝试了一下x^2+y^2==(Arcsin[Sin[x]])^2,积分结果仍然是0 即使Region输出的结果貌似是完全相同的 Clear["Global`*"] surfacea = ImplicitRegion[ x^2 + y^2 == (Cos[z])^2 && 0 < z < \[Pi]/2, {x, y, z}]; surfaceb = ImplicitRegion[ x^2 + y^2 == (ArcSin[Sin[z]])^2 && 0 < z < \[Pi]/2, {x, y, z}]; surfacec = ImplicitRegion[x^2 +

SetAttributes[{a, b, q}, Constant]; eqn = {Laplacian[w[x, y], {x, y}] == v[x, y], Laplacian[v[x, y], {x, y}] == 7}; bcs = {w[0, y] == 0, w[x, 20] == 0, w[10, y] == 0, w[x, 0] == 0, Derivative[1, 0][w][0, y] == 0, v[x, 20] == 0, Derivative[1, 0][w][10, y] == 0, v[x, 0] == 0}; ufun = NDSolveValue[{eqn, bcs}, w, {x, 0, 10}, {y, 0, 20}] NDSolveValue::fembdnl: The dependent variable in (w^(1,0))[0,y]==0 in the boundary condition DirichletCondition[(w^(1,0))[0,y]==0,x==0.] needs to be linear. 怎么搞成线性的啊

求助，各位大佬，程序中间有没有错误，现在的情况是没有办法出图 代码如下 中间要对特征值进行处理，最后画图，原本的问题是特征值排序的问题，我用了一位大佬解答的程序片段，现在的情况是没有办法出图，不清楚是哪里出了问题，结果是后面的图片 p = K0^2 \[Omega] - Ke^2 \[Omega] + Kf^2 \[Omega] - Kv^2 \[Omega] + \[Epsilon] \[Omega]^2 - \[Omega]^3 - 2 Ke Kf Subscript[l, B] + 2 K0 Kv Subscript[l, B] + \[Epsilon] \!\(\*SubsuperscriptBox[\(l\), \(B\), \(2\)]\) - \[Omega] \!\(\*SubsuperscriptBox[\(l\), \

各位大佬，我在计算特征值时，出现错序 跳变问题，直接画图就能出现。两个都出现了，应该是同一个问题，劳烦各位大佬解答，代码不长，只是占的页数比较大，请见谅。 程序一 \[Gamma] = 1; Subscript[\[Gamma], c] = 0*\[Gamma]; Subscript[\[Gamma], d] = 2*\[Gamma]; Subscript[\[Gamma], a] = 0*\[Gamma]; Subscript[\[Gamma], b] = 2*\[Gamma]; k = 1; J = 1; \[CapitalDelta] = 0.01; \[CapitalPhi] = 2*\[Pi]; n = 0;(*AW=0;*) U = N[BesselJ[n, AW]]; W = N[BesselJ[n, AW]]; TE = 20; e1 = Collect[Eigenvalues[1/I ({ {-I*Subscript[\[Gamma], a]/2, -k

想请问一下吧友我想用mathematica的dsolve解微分方程组，我带入的参数（如Y（T））都是矩阵形式的，但是dsolve一直不输出，想请问下如何解决，是不是必须一条一条把方程列出来啊。

无论是内置的Eigenvalues，还是使用SchurDecomposition，总会莫名的改变原矩阵特征值的顺序。比如一个4阶方阵，test ={{600, 20, 20, 20}, {20, 400, 20, 20}, {20, 20, 800, 20}, {20, 20, 20,1000}}，对角原远大于非对角元，显然特征值的顺序应该是600左右的一个，400左右的一个，800左右的一个，1000左右的一个。可是Eigenvalues一定会从大到小牌序，用Schur分解的方法就更离谱了： test 1=SchurDecomposition[N@test][[1]]; Transpose@test1.test.test1 给出的顺序是400,600,1000,800，并且还会随着非对角元

我用mathematica13.0计算一个很简单的极限并要求给出step by step solution, 发现答案没错但给的过程有明显错误（划线处），怎么回事？ Limit[(2*Sin[p*Pi] - 2*p*Pi*Cos[p*Pi])/p^2, p -> Plus[Infinity]]

求！ a-y/x-y 是关于θ的值，要画的也是θ=a-y/x-y 的三维平面 <math xmlns='http://www.w3.org/1998/Math/MathML' mathematica:form='StandardForm' xmlns:mathematica='http://www.wolfram.com/XML/'> <mrow> <mi>Show</mi> <mo>[</mo> <mrow> <mrow> <mi>Plot3D</mi> <mo>[</mo> <mrow> <mrow> <mi>NF</mi> <mo>/.</mo> <mi>vall11</mi> </mrow> <mo>,</mo> <mrow> <mo>{</mo> <mrow> <mi>&#952;</mi> <

\[CapitalOmega] = ImplicitRegion[True, {{z, 0, 200}}]; \[CapitalGamma]1 = NeumannValue[0 , z == 0 || z == 200]; TraditionalForm[\[CapitalGamma]1] {nO, pO} = NDSolveValue[ { D[n[t, z], t] == D[n[t, z], {z, 2}] + \[Tau]/ nb \[Alpha] (1 - R) Pdensity/(Sqrt[2 \[Pi] ] \[Tau]L hv ) gauss[t, 0, \[Tau]L/\[Tau]] Exp[-\[Alpha] Ldi z] + Ldi e^2 nb Ldi/(k T \[Epsilon] \[Epsilon]b) (0 + eF0O[z]) D[ n[t, z] , z] + Ldi e^2 nb Ldi/(k T \[Epsilon] \[Epsilon]b) (p0[z] - n0[z] + 1 - E^(-2 u) + p[t, z] - n[t, z]) n[t, z] + \[CapitalGamma]1, D[p[t, z], t] == 0.1 D[p[t, z], {z, 2}] + \[Tau]/ nb \[Alpha] (1 - R) Pde

\[Phi]1 = D[\[Phi][x, y, t], x, x, t] \[Phi]2 = D[\[Phi][x, y, t], x, x, x] \[Phi]3 = D[\[Phi][x, y, t], x, x, y] \[Phi]4 = D[\[Phi][x, y, t], y, y, t] \[Phi]5 = D[\[Phi][x, y, t], y, y, x] \[Phi]6 = D[\[Phi][x, y, t], y, y, y] \[Phi]7 = D[\[Phi][x, y, t], t] \[Phi]8 = D[\[Phi][x, y, t], x] \[Phi]9 = D[\[Phi][x, y, t], y] \[Phi] = \[Psi]0[y, t] + \[Psi] \[Phi]1 + \[Phi]4 - \[Lambda]^2 \[Phi]7 - \[Phi]9*\[Phi]2 + \[Phi]9*\ \[Phi]5 + \[Lambda]^2*\[Phi]9*\[Phi]8 + \[Phi]8*\[Phi]3 + \[Phi]8*\ \[Phi]6 - \[Lambda]^2*\[Phi]8*\[Phi]9 + \[Beta]*\[Phi]8 == 0 TraditionalForm[\[Phi]1 + \[Phi]4 - \[Lambda]

\[lbk]Omega[rbk]m = \[lbk]Pi[rbk]*4.3*10^5; m = 1.45*10^-7; k1 = k2 = \[lbk]Omega[rbk]m; \[lbk]CapitalDelta[rbk] = 1.2*\[lbk]Omega[rbk]m; \[lbk]CapitalDelta[rbk]2 = 2*\[lbk]Omega[rbk]m; P = 2*10^5; e = 1.6*10^-19; J = 1.5*k1; \[lbk]CurlyTheta[rbk] = 0; G = 0; \[lbk]HBar[rbk] = 1.05*10^-34; L = 0.025; \[lbk]Omega[rbk]L = (6*\[lbk]Pi[rbk]*10^8)/(1.064*10^-6); g = (\[lbk]CapitalDelta[rbk] + \[lbk]Omega[rbk]L)/(L - Q); Q = (\[lbk]HBar[rbk]*g*Abs[lbk]A[rbk]^2)/(m*\[lbk]Omega[rbk]m^2); \[lbk]Omega[rbk]1 = ((\[lbk]CapitalDelta[rbk] + \[lbk]Omega[rbk]L)*L)/(L - Q); \[lbk]CurlyEpsilon[rbk] = Sqrt[lbk](

如题。

mma 同样作为 lisp-like 语言有 Throw/Catch 有 Dialog 甚至有 Goto 却没有 call/cc 的官方实现也是够奇怪的

在尝试求解把vk/vkdig这样的变上限积分携带自变量代入到微分方程中求解时，发现它的计算速度非常的慢。请问该怎么优化呢？

Clear["Global`*"] Clear@Derivative DSolve[-g m + k 2 x[t] - 2 k l x[t]/Sqrt[l^2 + x[t]^2] + m (x^\[Prime]\[Prime])[t] == F[t], x[t], t] 求解微分方程的代码，但是运行后却直接把原式还给我了？不知道怎么回事？求大佬教教

链接不发了。楼下简述更新内容

最近在重看有限元相关的东西，结果在看到等参元这部分内容的时候发现了一个数学上的小问题，姑且拿出来水一帖。考虑下面这个参数方程： base = DeleteCases[{-(1/4) (-1 + s) (-1 + t) (1 + s + t), -(1/ 4) (1 + s) (-1 + s - t) (-1 + t), 1/4 (1 + s) (1 + t) (-1 + s + t), 1/4 (-1 + s) (1 + s - t) (1 + t), 1/2 (-1 + s^2) (-1 + t), -(1/2) (1 + s) (-1 + t^2), -(1/2) (-1 + s^2) (1 + t), 1/2 (-1 + s) (-1 + t^2)} /. t -> -1, 0][[{1, 3, 2}]] (* {1/2 (-1 + s) s, 1 - s^2, 1/2 s (1 + s)} *) 这个方程有如下性质：s为-1时，第1项为1，其余

其实本来想用的标题是《【交流】“这些问题你换软件也没用！”——受限于数学发展而无法或很难求解却经常有人想用Mathematica解的问题》。Mathematica作为一个软件，必然无法超越现有的数学理论，也就是说，对于尚未在理论上得到很好解决的问题，Mathematica是不太可能直接求解出来的——这话说出来大家都明白，但是，一般用户常常很难分辨哪些问题目前在理论上是很难或无法求解的，故开此交流帖，大家可以来谈谈自己所知道的相关问题。为避免

比如Reduce[BesselK[1, a]^2 == b && a > 0 && b > 0, a, Reals] 会提示无法利用 Reduce 现有的方法求解该系统

考虑到此类问题和“Mathematica已经装好了，但不能用”即“Mathematica一开始就不能用”有一定区别，特开专帖讨论。今后此类问题原则上请不要单开主题。 这里的“不能用”，包括不限于 1. Mathematica突然无法正常打开 2. Mathematica突然变得无法正常工作（1+1也算不了等） 3. Mathematica突然在执行正确的代码时出现莫名其妙的警告。 如果你的“Mathematica一开始就不能用”，请看置顶。 闲话少叙，说一说可能的修复方法： 此类问题似乎和Mathematica的Paclet更新不

c = 0.5 \[Alpha] = 0.9 d = 4 k = 4 \[Beta] = (d*(4 - \[Alpha])*(\[Alpha] - 2) + 4*(3 w + c*(\[Alpha] - k)) - \[Alpha]*(c*k + w - \[Alpha]))/(c* k*(4 - 3*\[Alpha])) t = (3*\[Alpha] + w*(8 - \[Alpha]) + d*(4 - \[Alpha]) - c*(4*k - 4 - \[Alpha]) - 4)/(4 - 3*\[Alpha]) Plot[{(2*w + \[Alpha]*(1 + c) + 2*k*c*(\[Beta] - 1))/(\[Alpha]*(4 - \[Alpha])), (2* w + \[Alpha] + \[Alpha]*c + \[Alpha]*t - 2*c*k)/(\[Alpha]*(4 - \[Alpha]))}, {w, 3.7, 4.8}] 这段代码有啥错误吗，我的输出（白色底）和原文（绿色底）不一样

大佬帮忙看看 Simplify[(((I Subscript[C, 1])/( 1 - x Subscript[C, 1] Subscript[L, 1]) + (I Subscript[C, 2])/( 1 - x Subscript[C, 2] Subscript[L, 2]) + (I Subscript[C, 3])/( 1 - x Subscript[C, 3] Subscript[L, 3]))*(I Subscript[C, 4])/( 1 - x Subscript[C, 4] Subscript[L, 4]))/(((I Subscript[C, 1])/( 1 - x Subscript[C, 1] Subscript[L, 1]) + (I Subscript[C, 2])/( 1 - x Subscript[C, 2] Subscript[L, 2]) + (I Subscript[C, 3])/( 1 - x Subscript[C, 3] Subscript[L, 3])) + (I Subscript[C, 4])/( 1 - x Subscript[C, 4] Subscript[L, 4])) + I x Subscript[C, 5], x] 怎么将这个公式化简通分成分

新人不知道该如何在最后一个输出里保留Pi。。。多谢了！

按照官方教程先pip install zmq了(图1)，然后 FindExternalEvaluators[lbk]"Python"[rbk] 没有结果， 又 FindExternalEvaluators[lbk][rbk] 结果有一个Jupyter，但是显示False。(图二) 然后随便试了用Python语言，报错连接失败。(图三) 附: 使用的是Anaconda，平时用VS Code或Jupyter notebook都能正常运行Python代码。

为什么设定横坐标刻度后上面的刻度就没了 \[CapitalDelta] = 1.66; s = 1; s1 = -1; \[Gamma] = 75*10^(-3); h = 6.58*^-16; p = 8.617*10^-5; \[Mu] = 0.8397209424904093`; \[Mu]1 = 0.8388894343099993`; T = 300; e = Sqrt[1.44 10^-7]; Plot[Re[1/ 16 (1/(1 + E^((-\[Mu] + 1/2 (s \[Gamma] - 2000000000000 h \[Pi] \[Nu]))/( p T))) - 1/( 1 + E^((-\[Mu] + 1/2 (s \[Gamma] + 2000000000000 h \[Pi] \[Nu]))/( p T)))) (1 + (-s \[Gamma] + \[CapitalDelta])^2/( 4000000000000000000000000 h^2 \[Pi]^2 \[Nu]^2)) + 1/16 (1/( 1 + E^((-\[Mu]1 + 1/2 (s1 \[Gamma] - 2000000000000 h \[Pi] \[Nu]))/(p T))) - 1/( 1

如题，想要实现“矩阵直乘再乘积等于对应乘积的直乘”运算，可以通过定义 A_\[CircleTimes]B_ := KroneckerProduct[A, B]; kroneckerProductRule = (a_\[CircleTimes]b_) . (c_\[CircleTimes]d_) :> \ (a . c)\[CircleTimes](b . d); 实现。现在想将这个运算推广到多个矩阵，也即 (A \[CircleTimes] B \[CircleTimes]C) . (D \[CircleTimes] E \[CircleTimes] F) = （A . D) \[CircleTimes] (B . E) \[CircleTimes] (C . F)，尝试模式匹配后发现不成功，甚至两项矩阵都未能实现想要的结果。请问有没有大佬能解答一下。代码如下：

代码如下：FullSimplify[ ImplicitD[x - y, (x - y)/(2 z) == n[x]/p[k[x], z] - n[y]/q[k[y], z], x - y, z]] 想要求d(x-y)/dz，其中p和q是k和z的函数，n和k是x或y的函数，求大佬帮忙看看是不是有什么问题，谢谢

\[CapitalOmega] = ImplicitRegion[True, {{z, 0, 200}}]; \[CapitalGamma]2 = NeumannValue[0., x == 0 || x == 200] s = NDSolveValue[{D[n[t, z], t] - D[n[t, z], {z, 2}] - \[Tau]/ nb \[Alpha] (1 - R) Pdensity/(Sqrt[2 \[Pi] ] \[Tau]L hv ) gauss[t, 0, \[Tau]L/\[Tau]] Exp[-\[Alpha] Ldi z] == -Ldi e^2 nb Ldi/(k \ T \[Epsilon] \[Epsilon]b) (0 + eF0O[z]) D[n[t, z] , z] - Ldi e^2 nb Ldi/(k T \[Epsilon] \[Epsilon]b) (p0[z] - n0[z] 0 + 1 - E^(-2 u)) n[t, z] + \[CapitalGamma]2, n[-10 \[Tau]L/\[Tau], z] == 0 }, n, {t, -10 \[Tau]L/\[Tau], 5}, {z} \[Element] \[CapitalOmega], PrecisionGoal -> 35, Method ->

重积分我想用那个boole函数表示结果很多简单的都积不出来，第二形曲面积分完全不知道该怎么用mma来算

求助有关ListVectorPlot的一些问题，下面的123分别对应着三幅图产生的问题： 1.首先Rainbow最亮的颜色是红色，下面这串代码将大于0.5的部分全部设置为最亮的颜色红色，但其实对于EVPA中的最大值是3.6，这样导致形成的图像大面积的全是红色（存在的问题可能是数据中存在较大的值，而大部分值相对较小，那么图例的颜色范围和实际图像中的颜色会不一致，个人猜测，不一定。）我想要的结果是只在最大值附近才是红色，让红色面积减少。 2.我通过将Vecto

现在有如图所示的数据点，z方向每隔0.05有一组数据分立分布的，每一组都形成了一条二维封闭曲线（原先数据量z方向上范围比较大，数据点已经够密了，为了减少贴出来的数据量截取了一小部分） 通过对数据的辐角和z坐标值排序，可以用ListLinePlot3D画出连成线的版本，效果还不错（部分地方变化较快有折线，可能也是后面画曲面有问题的原因） 想请问一下怎么让这些点（或二维线）生成一个三维面，类似一个筒的效果？ 用ListSurfacePlot3D画出来是很

各位大佬，求问

y = Integrate[Exp[x/2 (1 + x^2 - 2 x)], x]; Table[{x, y}, {x, 0, 1, 0.2}] 如例子所示，想求这个函数y(x)的数值Table，但是y是x的一个积分形式没法解析求解，我想生成一个{x,y}的table画它的函数图像但是用Integrate的话它不能直接取值会变成图里这个样子，求问一下怎样才能让它积分后再取值啊用NIntegrate的话它又不支持不定积分一定要输入上下限，如果能直接y=NIntegrate[y(x),x]就好了

如何使得f[g[x_]]的定义始终优先于g[x] 如上，我期望f[g[2]]始终输出4 代码如下 ClearAll["Global`*"] f[g[x_]]^:=2x; f[g[2]] g[x_]:=x f[g[2]]

现在是手动二分法，输入两个自变量x,y，寻找输出的两个量为特定的值时与之对应的xy，现在想能不能在mma里或者python中二分法求解x,y。已知函数单调，和xy大致范围，不会出现不同xy对应同一个值

如题，矩阵mA，mC是两个明确定义的矩阵；mB，mD是未定义的符号。请问如何让mA.mC \[CircleTimes] mB.mD在输出时让mA.mC为矩阵格式而非大括号格式？尝试了MatrixForm命令，但未能达成目的。将mA和mC矩阵化导致不计算矩阵乘积，而在最后使用MatrixForm命令没有结果。代码如下： Clear["Global`*"]; mA = {{1, 2}, {3, 4}}; mAform = mA // MatrixForm; mC = {{1, 0}, {0, 1}}; mCform = mC // MatrixForm; (mA\[CircleTimes]mB) . (mC\[CircleTimes]mD) /. (x_\[CircleTimes] y_) . (a_\[CircleTimes]b_) -> (x . a)\[CircleTimes]

下图是运算结果，我希望更改频率和线径来得出一个量来进行比较，假如如下图输入所示fs为3＊1的矩阵，线径为1＊3的矩阵，最后不应该α是一个3*3的矩阵，为啥是现在这种情况，请各位大佬帮忙看看 代码如下 fs = {70 Quantity[lbk]1, "Kilohertz"[rbk], 140 Quantity[lbk]1, "Kilohertz"[rbk], 150 Quantity[lbk]1, "Kilohertz"[rbk]};(*开关频率*) WireDs = {{0.1 Quantity[lbk]1, "Millimeters"[rbk], 0.07 Quantity[lbk]1, "Millimeters"[rbk], 0.05 Quantity[lbk]1, "Millimeters"[

后半截（不许问这个是什么或是怎么用）： s/1ytBfx9aZZLQvhSBvt9YqgA?pwd=wisz 两本书分别是 《Fundamental Finite Element Analysis and Applications: with Mathematica and Matlab Computations》 《Advanced Topics in Finite Element Analysis of Structures: With Mathematica and MATLAB Computations》 书的优点是讲解极为详细，几乎看不到跳步，即便撇开编程的部分也是很不错的有限元教材（尤其是第一册）。此外随书附赠的源代码（上面的链接中已包含）也是不错的学习材料。

c = 0.5 \[Alpha] = 0.9 d = 4 k = 4 \[Beta] = ( d (4 - \[Alpha]) (\[Alpha] - 2) + 4 [3 w + c (\[Alpha] - k)] - \[Alpha] (c k + w - \[Alpha]))/( c k (4 - 3 \[Alpha])) Plot[(2 w + \[Alpha] (1 + c) + 2 k c (\[Beta] - 1))/(\[Alpha] (4 - \[Alpha])), {w, 3.7, 4.8}]