Advanced Search

Question: Singularity

Question:Singularity

Stephen Moynihan 48 Maple
ode dsolve numeric differential-equations
March 12 2007
0
Hi, I have a ‘slight’ problem (you will probably recognize it Joe! :-) ). It concerns the values of Tau, omega in my worksheet (see below). If I set Tau=0.7, omega=0.7*m*Eta everything is rosy, and works fine. If I start tweaking these values (which I have to) things go a bit pair shaped. I either get an error message after the first call to dsolve (e.g. when tau=0.5, omega=0.7*m*Eta) : "Error, (in dsolve/numeric/checksing) ode system has a removable singularity at r=1. Initial data is restricted to {Phi(r) = .20650095602297*diff(Phi(r),r)+.82088920025557e-1*I*diff(Phi(r),r)}" Or I get an error message after BC2 (e.g. when tau=0.6, omega=0.2*m*Eta) : "Error, (in Phi_z) cannot evaluate the solution further right of 1.0000000, probably a singularity" Now I know there is a singularity (these equations describe a Black Hole) but I should be able to solve these equations for tau = 0.5,0.7,0.9,0.99, and where the factor in omega can be between 0 and (less than) 1. I also have a bit of a problem with an fprint command in one of my try ‘loops’, where I try to print out r from 2..3 and the corresponding Zp values. This problem is probably down to my lack of experience of using fprint in Maple, rather than the presence of any singularities! :-) Have any ideas as to how I could work around these problems? (my worksheet, in the form of a procedure that is read into Maple, is below) Cheers, the help is much appreciated. #The call to read in, and execute, the procedures in Maple: > restart; > read "Z:\\BH_v2.txt"; > get1(); > get2(); #The txt file BH_v2.txt: with(plots): #define variables: N:=2: mu:=0: m:=4: l:=2: sigma:=1: Tau:=0.7: M:=1/(2*(1-(Eta)^2)): a:=Eta: Eta:=Tau*sqrt(N/(N+1)): omega:=0.7*m*Eta: g:=r->sqrt((1-2*M/r^(2*N)+2*M*a^2/r^(2*N+2))^(-1)): h:=r->sqrt(r^2*(1+2*M*a^2/r^(2*N+2))): Omega:=r->2*M*a/(r^(2*N+2)+2*M*a^2): f:=r->r/(g(r)*h(r)): #define the Potential: V0:=r->(f(r))^2*sqrt(h(r))/(r^(N+1))*diff((((f(r))^2)*h(r)/r)*diff(sqrt(h(r))*r^N,r),r): V:=r->V0(r)+(f(r))^2*mu^2-(omega-m*Eta)^2+((f(r))^2/(r^2))*(l*(l+2*N)-m^2*(1-(r^2/(h(r)^2)))+4*(1-sigma)*((h(r))^2/(r^2)-1)): #define the ODE: ODE := -f(r)*(g(r)^(-1))*diff(f(r)*(g(r)^(-1))*diff(Phi(r),r),r)+2*I*(omega-m*Eta)*f(r)*(g(r)^(-1))*diff(Phi(r),r)+Phi(r)*(V(r)+(omega-m*Eta)^2)=0: #Define Boundary Conditions: YY:=diff(V(r),r): XX:=eval(YY,r=1): alpha1:=diff(g(r)^(-2),r): alpha:=alpha1/h(r): x:=r->int(g(r)/f(r),r): FF:=(XX*(r-1))/(alpha^2-2*I*alpha*(omega-m*Eta)): GG := subs(r=1.000005,FF): BC := ( Phi(1)=1, D(Phi)(1) = GG ): #procedure - get1 - calculates numerical solution of ODE # and then exports results to an external file: get1:=proc() #Declare variables: global dsol, dsol_1, file, tt, fd, Phi_z, DPhi_z: #Find numerical solution of system of differential equations: dsol := dsolve({ODE,BC},numeric, range=1..2); dsol_1 := dsolve({ODE,BC},numeric, output=listprocedure); #Open file to print numerical answer to: file:="Z:\\Black Holes\\Maple\\DATA\\BH_data.txt": try fd := fopen(file,'WRITE','TEXT'); for tt from 1 to 2 by 0.01 do fprintf(fd,"%a %a %a %a %a\n",eval([r,Re(Phi(r)),Im(Phi(r)),Re(diff(Phi(r),r)),Im(diff(Phi(r),r))],dsol(tt))[]) end do: finally close(file) end try: #define Phi(r) and DPhi(r) for use in solving Radial Eq Phi_z:=eval(Phi(r),dsol_1): DPhi_z:=eval(diff(Phi(r),r),dsol_1): #end procedure: end proc: #define second ODE: ODE2 := -(f(r)/g(r))*diff((f(r)/g(r))*diff(Psi1(r),r),r) + V(r)*Psi1(r): #define Boundary Conditions: BC2 := ( Psi1(2) = Phi_z(2), D(Psi1)(2) = DPhi_z(2)-I*(omega-m*Eta)*(g(2)/f(2))*Phi_z(2) ): #procedure - get2 - calculates numerical solution of 2nd ODE # and then exports results to an external file: get2:=proc() #Declare variables: global dsol2, dsol_2, file, file2, ttt, tttt, fd, Psi_z, DPsi_z, Zp: #Find numerical solution of system of differential equations: dsol2 := dsolve({ODE2,BC2},numeric, range=2..10); dsol_2 := dsolve({ODE2,BC2},numeric, output=listprocedure); Psi_z:=eval(Psi1(r),dsol_2): DPsi_z:=eval(diff(Psi1(r),r),dsol_2): Zp:=r->abs( (DPsi_z(r)+I*omega*Psi_z(r)) / (DPsi_z(r)-I*omega*Psi_z(r)) ); #Open file to print numerical answer to: #file:="Z:\\Black Holes\\Maple\\DATA\\Zp_data.txt": #try – HAVE PROBS’ WITH THIS BIT # fd := fopen(file,'WRITE','TEXT'); # for tttt from 2 to 3 by 0.0005 # do # fprintf(fd,"%a %a\n",eval([r,Zp(r)],Zp(tttt))[]) # end do: #finally close(file) #end try: #Open file to print numerical answer to: file:="Z:\\Black Holes\\Maple\\DATA\\BH_data2.txt": try fd := fopen(file2,'WRITE','TEXT'); for ttt from 2 to 10 by 0.01 do fprintf(fd,"%a %a %a %a %a\n",eval([r,Re(Psi1(r)),Im(Psi1(r)),Re(diff(Psi1(r),r)),Im(diff(Psi1(r),r))],dsol2(ttt))[]) end do: finally close(file2) end try: #end procedure: end proc: #N.B. get1() has to run for get2() to work.
Please Wait...