@rachie90 The situation which is described in the paragraph following the one containing equation (2.2) in Tsuchiya's paper is exactly what is happening in your case: G'/G'' is going to 0/0 (where I'm abusing the derivative notation a bit, but you know what I mean). G'' is going to 0 a little faster, so this is manifested in the floating-point arithmetic as the values getting extremely large.

So, I guess that it's time to try equation (2.3). I think (although I'm far from sure) that this is equivalent to handling r and z together, so we're just dealing with gradients and Hessians of functions of two variables.

I guess that you saw that Tsuchiya does the catenoid problem later in the paper.

@rlopez  Thanks, it works. A plot shows clearly that -W is the correct derivative. Either

plot(-W(0,t), t= 0..1);

or

pds:-plot(-w(x,t), x=0, t= 0..1);

gives the plot of diff(u(x,t), x) for x=0.

With pdsolve, there seems to be a very subtle difference between what works in this context and what is "too far from a standard form".

@rlopez  Thanks, it works. A plot shows clearly that -W is the correct derivative. Either

plot(-W(0,t), t= 0..1);

or

pds:-plot(-w(x,t), x=0, t= 0..1);

gives the plot of diff(u(x,t), x) for x=0.

With pdsolve, there seems to be a very subtle difference between what works in this context and what is "too far from a standard form".

Thank you for your Answer. I was not familiar with the method of extracting the solution procedures. To answer the question that you pose rhetorically in the title of your Answer, there are two reasons:

1. It just seems more natural to me to work with the plots, because I am well acquainted with the internal structure of plots and rtables, but not at all well acquainted with the internal structure of the modules returned by pdsolve. Certainly this is just a matter of taste.
2. It hardly matters for this small example, but the plot method seems to use substantially less cputime (about 1/3 the time in this example). I haven't looked into why this is, but my guess is that the plots allow for more processing by evalhf in a "batch mode", i.e., fewer total calls to evalhf with each doing more work.

Robert Lopez wrote:

To get the values of the highest-ordered derivative, ... trick the system into thinking ... that the derivative is actually one of the system's dependent variables.

Could you give an example of this with pdsolve(..., numeric)? I am familiar with how to do this with dsolve(..., numeric), but when I try with pdsolve, I get an error message about the system being "too far from a standard form".

PDE:= {diff(u(x,t),t) = -diff(u(x,t),x), diff(u(x,t),x)= w(x,t)}:
IBC:= {u(x,0)=sin(2*Pi*x),u(0,t)=-sin(2*Pi*t)}:
pds:= pdsolve(PDE, IBC, numeric, time=t, range=0..1);
Error, (in pdsolve/numeric/par_hyp) input system is too far from a 'standard' form (see ?pdsolve,numeric for more detail)

Thank you for your Answer. I was not familiar with the method of extracting the solution procedures. To answer the question that you pose rhetorically in the title of your Answer, there are two reasons:

1. It just seems more natural to me to work with the plots, because I am well acquainted with the internal structure of plots and rtables, but not at all well acquainted with the internal structure of the modules returned by pdsolve. Certainly this is just a matter of taste.
2. It hardly matters for this small example, but the plot method seems to use substantially less cputime (about 1/3 the time in this example). I haven't looked into why this is, but my guess is that the plots allow for more processing by evalhf in a "batch mode", i.e., fewer total calls to evalhf with each doing more work.

Robert Lopez wrote:

To get the values of the highest-ordered derivative, ... trick the system into thinking ... that the derivative is actually one of the system's dependent variables.

Could you give an example of this with pdsolve(..., numeric)? I am familiar with how to do this with dsolve(..., numeric), but when I try with pdsolve, I get an error message about the system being "too far from a standard form".

PDE:= {diff(u(x,t),t) = -diff(u(x,t),x), diff(u(x,t),x)= w(x,t)}:
IBC:= {u(x,0)=sin(2*Pi*x),u(0,t)=-sin(2*Pi*t)}:
pds:= pdsolve(PDE, IBC, numeric, time=t, range=0..1);
Error, (in pdsolve/numeric/par_hyp) input system is too far from a 'standard' form (see ?pdsolve,numeric for more detail)

@J4James I can't get it. The obvious thing to try is

pds:-plot([diff(u(x,t), x)], x=0, t= 0..1, numpoints= 50);

But that just gives a constant 0 plot, which is obviously wrong.

@J4James I can't get it. The obvious thing to try is

pds:-plot([diff(u(x,t), x)], x=0, t= 0..1, numpoints= 50);

But that just gives a constant 0 plot, which is obviously wrong.

Do you have a reference on this "relaxation method" that I could look at? Possibly something on the web? I looked at the Wikipedia page "relaxation (iterative method)", but I couldn't relate it to your situation. How is the relaxation factor 1.2 chosen?

Your iteration uses r and z; why is theta not also used? Why is the iteration not done in Cartesian coordinates?

I don't think that it makes sense to sort the points. We used that code when we were dealing with a surface, a cylinder, of points. This time it is just a curve. They are already in a natural order.

You should evalf the initial points.

@kswong89 There are so many problems with this....

First, since you have F:= ... rather than F = ..., your expression is an assignment; it is not an equation. Since it is not an equation, it doesn't make sense to "solve" it. But I think that you mean "evaluate" rather than "solve".

Second, m and n are bound variables (as opposed to free variables), since they are index variables of summations. It doesn't make sense to assign values to them, or even to refer to them in any context more global than the summations, such as the overall expression.

Third, if you set the other variables to 1 that you wanted, two of the sines will be 0 for any m and any n. Thus, you'd be dividing by 0.

Fourth, I cannot expand your .png thumbnail to view it. This may be a limitation of the ability of my software. If it was .jpg, I could expand it. Anyway, please post full-size images.

@kswong89 There are so many problems with this....

First, since you have F:= ... rather than F = ..., your expression is an assignment; it is not an equation. Since it is not an equation, it doesn't make sense to "solve" it. But I think that you mean "evaluate" rather than "solve".

Second, m and n are bound variables (as opposed to free variables), since they are index variables of summations. It doesn't make sense to assign values to them, or even to refer to them in any context more global than the summations, such as the overall expression.

Third, if you set the other variables to 1 that you wanted, two of the sines will be 0 for any m and any n. Thus, you'd be dividing by 0.

Fourth, I cannot expand your .png thumbnail to view it. This may be a limitation of the ability of my software. If it was .jpg, I could expand it. Anyway, please post full-size images.

@Markiyan Hirnyk I see. Thank you.

@Markiyan Hirnyk I see. Thank you.

@Markiyan Hirnyk So isn't that worth a vote up?

@Markiyan Hirnyk So isn't that worth a vote up?

I think that the results of your isolve are difficult to interpret when compared to the corresponding results from msolve:

msolve(2275*s=2503, 1563);
                           {s = 976}

I don't see the significance of the chrem answer. Could you explain that?