求助:利用FindRoot求解自定义函数的零点
mathematica吧
全部回复
仅看楼主
level 7
Cauchy⚡ 楼主
我想写一个自定义函数(比如f(x)),然后用FindRoot来求这个函数的零点,但是我想不到如何实现我想要的效果。具体如下。
第一个问题。我想要让FindRoot在每一步迭代都print出目前的x和目前函数f(x)的值,这样我能清晰地看出寻根的进程究竟如何,不然的话我只能一直等待漫长的求解过程而且还不知道是不是求解遇到了什么困难。
我看到FindRoot有个选项EvaluationMonitor有类似的效果,但这个选项只能在方程求解完毕后才输出x的变化,所以应该满足不了我的需求。
同样的问题放在Matlab上,就很容易解决。只需要写一个函数文件,在这个函数文件内部加上display(x)和display(f(x))就可以了。但是我不清楚MMA如何实现这个功能。因为MMA给我的印象一直都是函数式编程,我不太清楚在MMA中如何写出灵活的程序实现我的需求。
第二个问题。我先贴上我的代码:
Clear["Global`*"];
Clear[Derivative];
$MaxPrecision=Infinity;
$MinPrecision=25;
$Assumptions=L>0&&L<1;
\[Alpha]=1`25;
J[L_]:={{0,-L-I Sqrt[1-L^2],0},{I/2,-\[Alpha],-(I/2)},{0,-L+I Sqrt[1-L^2],0}};
zStar[L_]:=I L-Sqrt[1-L^2];
zbarStar[L_]:=-I L-Sqrt[1-L^2];
T[L_]:=Module[{T0=Chop@Transpose[Eigensystem[J[L]][[2]]]},Dot[T0,{{1,0,0},{0,Exp[I x]/.Quiet@Solve[{Exp[I x]T0[[3,2]]==Exp[-I x]Conjugate[T0[[1,2]]]},x][[1]],0},{0,0,1}}]];
f[z_,y_,zbar_,L_]:={I z y,L-(z-zbar)/(2I)-\[Alpha] y,-I zbar y}
\[Alpha]1=(-0.0569346681735483100196233345904580474`25.+0.03380509290798852465990285780535192822`25. I);
\[Beta]1=1;
\[Alpha]2=1;
\[Beta]2=(-0.5621728742330930
18969289235
57469160676`25.-0.15431958634948985984144944464738474976`25. I);
\[Alpha]3=(1.12190799111714769562606230101308634098`25.-0.666
13550391341
96995510221227089331722`25. I);
\[Beta]3=0.7671179743387728812498376749441521484`25.;
u0[A_,B_]:=((\[Alpha]1 A^2+\[Beta]1 B)/2-(\[Alpha]1 A^2-\[Beta]1 B)/2 (2k-1))k^2 (1-k);
v0[A_,B_]:=((\[Alpha]2 A^2+\[Beta]2 B)/2-(\[Alpha]2 A^2-\[Beta]2 B)/2 (2k-1))k (1-k)^2;
w0[A_,B_]:=((\[Alpha]3 A^2+\[Beta]3 B)/2-(\[Alpha]3 A^2-\[Beta]3 B)/2 (2k-1))k^2 (1-k)^2;
dkdt=k(1-k);
g[u_,v_,w_,L_]:=Dot[Inverse[T[L]],f[({zStar[L],0,zbarStar[L]}+Dot[T[L],{u,v,w}])[[1]],({zStar[L],0,zbarStar[L]}+Dot[T[L],{u,v,w}])[[2]],({zStar[L],0,zbarStar[L]}+Dot[T[L],{u,v,w}])[[3]],L]]
targetN[L_,A_,B_]:={Chop@Expand@Integrate[Chop@Apart@(g[u0[A,B],v0[A,B],w0[A,B],L][[1]]/dkdt),{k,0,1}],Chop@Expand@Integrate[Chop@Apart@(g[u0[A,B],v0[A,B],w0[A,B],L][[2]]/dkdt),{k,0,1}],Chop@Expand@Integrate[Chop@Apart@(g[u0[A,B],v0[A,B],w0[A,B],L][[3]]/dkdt),{k,0,1}]};
FindRoot[{targetN[L,A,B]=={0,0,0},L>0&&L<1,A>0,B>0},{{L,0.9},{A,5},{B,5}}]
这段代码比较繁杂,不过关键的问题出在最后的
FindRoot[{targetN[L,A,B]=={0,0,0},L>0&&L<1,A>0,B>0},{{L,0.9},{A,5},{B,5}}]
这串代码。
大致就是我先定义了一些变量和函数,然后利用这些来定义我最终需要求解的函数targetN. L的具体取值会影响计算中一些运算的结果(比如EigenSystem),所以我不能写出targetN的显式表达式(而且显式表达式也极端复杂),只好每给定一个L(以及A, B),计算targetN[L,A,B]的值,以此来求targetN的零点。但是如果我运行
FindRoot[{targetN[L,A,B]=={0,0,0},L>0&&L<1,A>0,B>0},{{L,0.9},{A,5},{B,5}}]
的话,会报错:
Integral of (b^2 <<4>>)/((1.0-1.000000000000000000000000 k)^2 (-1.000000000000000000000000+L^2)) does not converge on {0,1}.
也就是说targetN[L,A,B]先把L, A, B当作没有赋予具体数值的符号来计算targetN的值,然后报错说积分不收敛,这恰恰是我想要避免的。我想要的是每代入一个具体的L, A, B的值都能给出相应的targetN[L,A,B]的值,利用牛顿法一步步迭代下去,这样就能避免符号计算带来的麻烦,至少Matlab中是这么一个过程。
所以第二个问题大概就是这个。
2020年01月02日 10点01分 1
level 7
Cauchy⚡ 楼主
报错中的b应该是B,在此纠正一下:
Integral of (B^2 <<4>>)/((1.0-1.000000000000000000000000 k)^2 (-1.000000000000000000000000+L^2)) does not converge on {0,1}.
2020年01月02日 10点01分 2
吧务
level 10
首先,EvaluationMonitor是可以算一个输出一个的,你可能没用:>而用了->导致包含Print的表达式预先被计算为不包含Print的表达式。
然后我估计你报错是因为你输入的方程具有targetN[L, A, B] == {0, 0, 0}形式,系统必须对左端预先进行代数化简来确认是否可与右端由Equal来Thread,此时是不代入具体数值的,这时系统发现1/dkdt不被消去而又积分发散,就报错了。
2020年01月02日 11点01分 3
level 7
Cauchy⚡ 楼主
@asdasd1dsadsa EvaluationMonitor那个的确是我疏忽了。多谢!
至于第二个问题,我将targetN[L, A, B]拆成了三个函数
targetN1[L, A, B],targetN2[L, A, B],targetN3[L, A, B],
也就是targetN[L, A, B]的三个分量,分别让它们等于0,也就是
FindRoot[{targetN1[L,a,b]==0,targetN2[L,a,b]==0,targetN3[L,a,b]==0},{{L,0.6},{a,5},{b,5}},EvaluationMonitor:>Print["L = ",L," a = ",a," b = ",b," targetN =",targetN[L,a,b]]]
但是仍然报错,比如
Integral of (B^2 <<3>> ((0.*10^-3+7.42 I)-(0.*10^-3+<<43>> I) k+(<<1>>) <<1>>+(<<46>>-<<44>> I) (-<<1>>)^2.000000000000000000000000-((0.077
15979317474
492670935610-0.2810864371165465090641675 I) Sqrt[<<21>> L-<<1>>])/Sqrt[1.00000000000000000000 L+1.00000000000000000000 I Sqrt[Plus[<<2>>]]]))/((1.0-1.000000000000000000000000 k)^2 (-1.000000000000000000000000+L^2)) does not converge on {0,1}.
这说明L仍然被用于符号计算了。
然而,如果我将L当作一个常量L0=0.964327
18969383968
14391581`25(这个L0是有依据的,不是瞎取的),只利用targetN1和targetN2来解a和b,也就是
FindRoot[{targetN1[L0,A,B]==0,targetN2[L0,A,B]==0},{{A,5},{B,5}},EvaluationMonitor:>Print["A = ",A," B = ",B," targetN =",targetN1[L0,A,B]]]
则可以解出解
{A->2.1068 +0.180286 I,B->-0.0452825-0.0855974 I}
如果我将a和b都定为5,将L当作变量,只用targetN1[L,5,5]一个函数来求解L,却会报错说积分不收敛。
所以似乎L, A, B在函数targetN1的求解过程的地位并不“等价”?似乎只要L被包含在求解过程中,求解就会报错。
2020年01月02日 13点01分 4
我所知道/认为的是:报错应该是在代数化简的时候报的,但是能数值上解出来——FindRoot也确实返回了结果。
2020年01月03日 05点01分
@asdasd1dsadsa 如果是这样的话,我觉得可能是我定义的函数targetN中嵌套了很多其他函数,使得MMA先去化简这些表达式(比如积分),才导致了这个情况。MMA有没有办法像Matlab那样单独写一个m文件来当作函数调用呢?这样也许能避免过多嵌套。
2020年01月03日 07点01分
@asdasd1dsadsa 我查到的资料基本都是用自定义函数来定义一些简单的函数。
2020年01月03日 07点01分
@Cauchy⚡ 避免代数化简应该有别的办法,但我不知道怎么搞。
2020年01月03日 07点01分
吧务
level 15
1. 捧着个手机没法测代码所以我不敢说死,但我印象里编程层面上的计算次序问题是不会导致积分不收敛的,正常应该是报被积函数非数值。你要好好检查下被积函数。
2. 数值计算一般用NIntegrate而非Integrate。
2020年02月01日 04点02分 5
1