我想写一个自定义函数(比如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.56217287423309301896928923557469160676`25.-0.15431958634948985984144944464738474976`25. I);
\[Alpha]3=(1.12190799111714769562606230101308634098`25.-0.6661355039134196995510221227089331722`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中是这么一个过程。
所以第二个问题大概就是这个。
第一个问题。我想要让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.56217287423309301896928923557469160676`25.-0.15431958634948985984144944464738474976`25. I);
\[Alpha]3=(1.12190799111714769562606230101308634098`25.-0.6661355039134196995510221227089331722`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中是这么一个过程。
所以第二个问题大概就是这个。