with(DEtools): plotsetup(x11): test0:=proc(Nu,RR) local F,F1,F2,G1,G2; F:=-(-r-I*x)^2*(p2+I*q2)+(-r-I*x)*I*(p1+I*q1)+((-r-I*x)^2-Nu^2)*(p+I*q); F:=subs(r=RR,F); F:=expand(F); print(F); # return(F); F1:=coeff(F,I,0); F2:=coeff(F,I,1); # print(expand(F-F1-F2*I)); G1:=solve({F1=0,F2=0},{p2,q2}); return(G1); end: vv:=proc(Nu,RR,X) local Rule,P; P:=-RR-I*X; Rule:={x=X, p=Re(evalf(HankelH2(Nu,P))), p1=Re(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH2(Nu,z)),x)))), p2=Re(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH2(Nu,z)),x$2)))), q=Im(evalf(HankelH2(Nu,P))), q1=Im(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH2(Nu,z)),x)))), q2=Im(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH2(Nu,z)),x$2)))) }; return(Rule); end: ## test0 で導出した方程式が正しいか数値的に調べる. debug 方法1. test0a:=proc() local A,R,RR,Nu; Nu:=1/3; RR:=10; A:=test0(Nu,RR); R:=vv(Nu,RR,-2); print(A); print(R); return(subs(R,A)); end: ## 11.11 9:46 右=左 を戻すようになった. OK. ## test0a() がおかしいので, その原因を探る... vv0(1/3,2,-2) ## vv を修正. cf. vv0 in Buggy/han.ml-1110 test1:=proc() local Eqs,G1,Ans,P,Nu,RR,PP,Rule; Nu:=1/3; RR:=2; PP:=-2; P:=-RR-PP*I; Rule:=vv(Nu,RR,PP); G1:=test0(Nu,RR); G1:=subs(p=f0(x),p1=f1(x),q=g0(x),q1=g1(x),G1); print(G1); Eqs:={diff(f0(x),x)=f1(x),diff(f1(x),x)=p2, diff(g0(x),x)=g1(x),diff(g1(x),x)=q2, f0(PP)=subs(Rule,p),f1(PP)=subs(Rule,p1), g0(PP)=subs(Rule,q),g1(PP)=subs(Rule,q1) }; Eqs:=subs(G1,Eqs); print(Eqs); Ans:=dsolve(Eqs,numeric,range=-2..2); return(Ans); end: hc:=proc(Nu,z) local V; V:=evalf((1+exp(2*Pi*I*Nu))*HankelH1(Nu,z)-HankelH2(Nu,z)); return(V); end: # 漸近展開を使う. hc0:=proc(Nu,z) local V; V:=evalf((1+exp(2*Pi*I*Nu))*hc1(Nu,z,10)-hc2(Nu,z,10)); return(V); end: ## HankelH1 |Z|>=10, N=10 hc1:=proc(Nu,Z,N) local V,V0,i; V0:=evalf((2/(Pi*Z))^(1/2)*exp(I*(Z-Pi*Nu/2-Pi/4))); V:=1; for i from 1 to N do V:=V+evalf(pochhammer(1/2-Nu,i)*pochhammer(1/2+Nu,i)/ ((i!)*(2*I*Z)^i)); od: return (V0*V); end: ## HankelH2 |Z|>=10, N=10 hc2:=proc(Nu,Z,N) local V,V0,i; V0:=evalf((2/(Pi*Z))^(1/2)*exp(-I*(Z-Pi*Nu/2-Pi/4))); V:=1; for i from 1 to N do V:=V+evalf((-1)^i*pochhammer(1/2-Nu,i)*pochhammer(1/2+Nu,i)/ ((i!)*(2*I*Z)^i)); od: return (V0*V); end: vvh1:=proc(Nu,RR,X) local Rule,P; P:=-RR-I*X; Rule:={x=X, p=Re(evalf(HankelH1(Nu,P))), p1=Re(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH1(Nu,z)),x)))), p2=Re(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH1(Nu,z)),x$2)))), q=Im(evalf(HankelH1(Nu,P))), q1=Im(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH1(Nu,z)),x)))), q2=Im(evalf(subs(x=X,diff(subs(z=-RR-I*x,HankelH1(Nu,z)),x$2)))) }; return(Rule); end: # H1 の値を argz >= pi で計算してみる. A:=test1h(); A(2); どうやらこれでOK. test1h1:=proc() local Eqs,G1,Ans,P,Nu,RR,PP,Rule; Nu:=1/3; RR:=10; PP:=-2; P:=-RR-PP*I; Rule:=vvh1(Nu,RR,PP); G1:=test0(Nu,RR); G1:=subs(p=f0(x),p1=f1(x),q=g0(x),q1=g1(x),G1); print(G1); Eqs:={diff(f0(x),x)=f1(x),diff(f1(x),x)=p2, diff(g0(x),x)=g1(x),diff(g1(x),x)=q2, f0(PP)=subs(Rule,p),f1(PP)=subs(Rule,p1), g0(PP)=subs(Rule,q),g1(PP)=subs(Rule,q1) }; Eqs:=subs(G1,Eqs); print(Eqs); Ans:=dsolve(Eqs,numeric,range=-2..2); return(Ans); end: # これは不連続. arg z = pi に branch cut が入っている. # Try also plot(Im(sqrt(-2-x*I)),x=-2..2); # han-2hg-out.png h2g:=proc() plot(Im(HankelH2(1/3,-10-I*x)),x=-2..2); end: # これは連続. h2g0:=proc() plot(Im(HankelH2(1/3,2-I*x)),x=-2..2); end: test1junk:=proc() DEplot(-(1+I*x)^2*diff(y(x),x$2)-(1+I*x)*I*diff(y(x),x)+((1+I*x)^2-(1/3)^2)*y(x)=0, y(x), x=-1..1,[[y(0)=3.753,D(y)(0)=-1.40143]]); end: # See dsolve[numeric] # ss:=test0x() # ss(0); ss(3.14); test0x:=proc() RETURN( dsolve({diff(y0(x),x)=y1(x), diff(y1(x),x)=-y0(x),y0(0)=1, y1(0)=0}, numeric, range=0..5) ); end: