## 30 Reputation

2 years, 310 days

## Re: Caption?...

Thank you very much for your help. This code works well.

## A problem on Finding singular points...

Hi

Hope this message finds you well.

I have just found your code on finding the singular point of a differential equation.

However, I failed to generalize the code to the case in which the sol is expressed as a function of free parameters.

Please take a look at the file attached and I would appreciate for your help to overcome this problem.

My best regards

A_question_on_find_sing-08-02-2020.mw

 > restart;
 >
 > ode12:=(a,b,c0,c1)->    {-0.1*a*diff(u(z),z,z) +      (z-2*diff(v(z)^(-1/2),z))*diff(u(z),z)+(3-2*b*diff(v(z)^(-1/2),z,z))*u(z)= 0, 0.1*diff(v(z),z,z)+0.01*z*a*diff(v(z),z)+0.02*v(z)-u(z)*b*v(z)^(1/2)=0,u(0)=c0, v(0)=0.1, D(u)(0)=c1, D(v)(0)=0}; (1)
 > sol:=(a,b,c0,c1)-> dsolve(ode12(a,b,c0,c1), numeric):
 > sol(1,1,0,1); (2)
 > sol(1,1,0,1)(0.2); (3)
 > FindSingularity:= proc(sol(1,1,0,1), max::numeric)      try sol(1,1,0,1)(max)      catch "cannot evaluate the solution further":           return lastexception[-1]      end try;      FAIL end proc:
 > FindSingularity(sol, 4); (4)
 > FindSingularity(sol, 10); (5)
 > How can I fix the error?
 > #How can I define the FindSingularity procedure as a function of a,b,c0,c1? Then defining a loop via for cammand in order to get different singular points and their corresponding values for u(z) and (z), for different values of a,b,c0,c1?
 > #For exmple:
 > for a from 0.1 by 0.1 to 1 do for b from -1 by 0.1 to 1 do print(a, b, u(at singular z),v(at singular z)) od; od;
 > Would you be so kind as to let me know what does -1 means within the command: return lastexception[-1]
 > > ## Re: What do you hope to achieve??!!!...

Hi

Thank you very much for the time you devoted to this issue.

My problem is solved and I think that your code works well.

My best regards.

## > ....

 > restart;
 > # eq1 is the Non-linear differential equation and eq2 is the linear one that can be obtained by neglecting the quadratic terms in eq1.
 > eq1:=(a,b)->-(2/3)*b*(1+y(x))*(1+x)*(1+a*(1+x)^3)^2*(diff(y(x), x, x))+(8/9)*b*(1+a*(1+x)^3)^2*(diff(y(x), x))^2-(1/3)*b*(1+a*(1+x))*(1+y(x))*a*(1+a*(1+x)^3)*(diff(y(x), x))+y(x)*(1+y(x))^2*(1+x)^2*a^2*(2+a*(1+x)); (1)
 > eq2:=(a,b)->-(2/3)*b*(1+a*(1+x)^3)^2*(1+x)*(diff(y(x), x, x))-(1/3*(1+a*(1+x)))*a*(1+a*(1+x)^3)*b*(diff(y(x), x))+a^2*y(x)*(2+a*(1+x))*(1+x)^2; (2)
 > s1:=(a,b,y0,x0,x1)->dsolve({eq1(a,b),y(x0)=y0,(D(y))(x0)=-y0/(1+x0)},type=numeric, output = listprocedure, range=x0..x1); s2:=(a,b,y0,x0,x1)->dsolve({eq2(a,b),y(x0)=y0,(D(y))(x0)=-y0/(1+x0)},type=numeric, output = listprocedure, range=x0..x1);  (3)
 > s1(0.322,2.05e-05,0.01,5,-0.6)(4.0396832);
 > s1(0.322,2.05e-05,0.01,5,-0.6)(4.0396833); (4)
 > s2(0.322,2.05e-05,0.01,5,-0.6)(4.0396832); (5)
 > s11:=(a,b,y0,x0,x1,xx)->eval(y(x),s1(a,b,y0,x0,x1))(xx); s22:=(a,b,y0,x0,x1,xx)->eval(y(x),s2(a,b,y0,x0,x1))(xx);  (6)
 > # To obtain a certain value of y(x) from eq1 and eq2:
 > s11(0.322,2.05e-05,0.01,5,-0.6,4.0396833); s22(0.322,2.05e-05,0.01,5,-0.6,4.0396833);  (7)
 > # Now, how can I find those values of x at which the Non-linear solution is singular and diverges? #To this aim I tried to make a procedure:
 > y1:=proc(a,b,y0,x0,x1,xx) if s11(a,b,y0,x0,x1,xx)>=10^10 then print(s11(a,b,y0,x0,x1,xx), s22(a,b,x1,x0,x1,xx), y0, xx) end if; end proc; (8)
 > y1(0.322,2.05e-05,0.01,5,-0.6,4.0396833); (9)
 > # Instead of the condition, if s11(a,b,y0,x0,x1,xx)>=10^10, can we use the command singular (or other commands that I do not know) to find out at which value of x the Non-linear solution diverges? A command that checks wheather the solution is singular and then reports the corresponding value of x.
 > #Now I want to make a loop on the initial value y0 and xx with for command so that 1) by running xx, the loop finds the singular point (here at x=4.0396833) for each y0, 2) evaluates the value of Linear solution at x=4.0396833 and then write them in the set [xx,yLinear(xx)] in a text file.
 > for y0 from 0.01 by 0.001 to 0.02 do for xx from 4.0396839 by -0.0000001 to 1.0396833 do try print(y1(0.322,2.05e-05,y0,5,-0.6,xx)) end try; next; od; od;              > #It seems that only y0=0.01 is used but the loop over y0 has not been continued.
 > >
 >  Sorry for the problem! The file has been attached to the message.

## A sample code ...

Hi

Thank you for your response.

I have provided a sample that illustrates the situation and problems with which I faced.

My best regards

## @Kitonum    Thank you so muc...

Thank you so much for your kind help and answer.

## @Kitonum  Many thanks for your com...

Best

Thank you so much Thomas. But as I write the command with(Interpolation): the following error appears:

Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received Interpolation

How can I deal with this problem?

Thank you so much and best regards.