## 7776 Reputation

17 years, 363 days

## Need the worksheet...

@simplevn1967 I cannot diagnose the problem without seeing your worksheet.

As you reply to this message, you will see a large green upward pointing arrow in the toolbar.  Click on it to attach your worksheet to your reply.

## When is a solution not a solution...

@mmcdara Your calculation imposes the initial condition y(0)=0 but the result does not satisfy y(0)=0.  In what way is that a solution to the initial value problem?

## Deciphering...

@Carl Love That's cool, but perhaps too clever for my taste.  I prefer more verbose code that one can verify at a glance rather a overly terse code that requires effort to decipher.   Vote up, nonetheless.

## Comment...

@mmcdara I had no intention of criticizing your suggestion.  My remark was directed to the OP, who in response to my question of whether he knew how to write the equation of a tangent plane to a sphere, he responded with an appeal to Maple's geom3d package, while the simple answer is just (X − P) . (P − C) = 0, where C is the sphere's center, P is the point of tangency, and X is an arbitrary point on the tangent plane.

## Triangulation of the domain is the bottl...

@rlopez Yes, I recall having seen finite difference schemes applied to non-rectangular domains. I don't think that's worth the effort.

Writing a finite element scheme using 1st degree elements will take just a few dozens lines of Maple code and it will work for arbitrary-shaped domains. That's provided that we have a general proc for triangulating the domain. Maple's GraphTheory:-DelaunayTriangulation comes close to fulfilling that task but is not quite adequate.  For instance, it is not suitable for triangulating nonconvex domains.  Even on convex domains Delaunay triangulation produces poor results for finite element applications.

Jonathan Shewchuck's Triangle library (written in C) produces perfect triangulations of arbitrary two-dimensional domains bounded by polygons and possibly containing one or more polygonal holes.  I have made extensive use of that code in my programs (also in C).  I suspect that it would be possible to compile Shewchuck's C code into a Maple library but I haven't looked into that.  Once that's done, a slew of possibilities will open for implementing finite element code in Maple.  Perhaps someone at Maplesoft would be interested in looking into that.

## Can you do it for a single point?...

It's not clear how much help you will need for this.  Do you know how to write the equation of the tangent plane to a sphere at a given point?  If you explain how you would do that, then I will be happy to show how to do that for a list of points.

## Rectangular domains...

@rlopez That's a nice video.  In it you explain a variety of ways, and their pros and cons, of solving elliptic and parabolic PDEs.

The finite differences method that lies at the foundation of all those methods, limits the spatial domain to a rectangle (of any dimension).  It will be difficult, if not impossible, to extend those methods to more general domains. (I am thinking of a square domain from which one or more holes have been punched out).  That's not a fault of your presentation; that's an inherent limitation of the finite difference methods.

At one point you say: "I have a feeling that this is why Maple does not have an elliptic numeric solver cause it is hard to have it robust enough to handle a variety of problems."

I agree.  Implementing a solver that works on arbitrary-shaped domains is a monumental task. There are specialized companies—Comsol and LS-Dyna come to mind—devoted to that task.  Their underlying method is finite elements rather than finite differences.  Perhaps one day Maple will branch out in that direction but I don't envision that happening anytime soon.

## @sand15 If the predator starts at a...

@sand15 If the predator starts at an arbitrary position, of course there is no reason for it to capture the target exactly at the origin.

The point of the puzzle is to find out where the predator is to start so that the capture occurs at the origin.  The answer is:  if the predator starts anywhere on an ellipse (that is to be determined), then the capture does occur at the origin.

## @sand15 Yes, you are quite correct....

@sand15 Yes, you are quite correct.  I had not explored the possibility of such small values of E1 and E2. Thanks for pointing that out.

## @raj2018 I am getting zero for f1 e...

@raj2018 I am getting zero for f1 evaluated at phi=0, as you expect.  Have a look.

 > restart;
 > f1:=(1-delta[h]-delta[b])*(1-phi/(kc-3/2))^(-kc+3/2)*(-kc+3/2)/((kc-3/2)*(1-phi/(kc-3/2)))+delta[h]*(1-phi/(sigma[h]*(kh-3/2)))^(-kh+3/2)*(-kh+3/2)/((kh-3/2)*(1-phi/(sigma[h]*(kh-3/2))))-(1/18)*delta[b]*(3*sqrt((M-u0b+sqrt(3)*sqrt(mu*sigma[b]))^2+2*mu*phi)*mu-3*sqrt((M-u0b-sqrt(3)*sqrt(mu*sigma[b]))^2+2*mu*phi)*mu)*sqrt(3)/(mu*sqrt(mu*sigma[b]))-(1/18)*(-3*sqrt((M+sqrt(3)*sqrt(sigma[i]))^2-2*phi)+3*sqrt((M-sqrt(3)*sqrt(sigma[i]))^2-2*phi))*sqrt(3)/sqrt(sigma[i]);

 > g1 := subs(phi=0, f1);

 > simplify(g1) assuming positive,         M - u0b - sqrt(3)*sqrt(mu)*sqrt(sigma[b]) > 0,         M - u0b + sqrt(3)*sqrt(mu)*sqrt(sigma[b]) > 0,         M - sqrt(3)*sqrt(sigma[i]) >0;

 >

## Where is the derivative?...

You say that "the eval command does not gives the correct value of its first derivative" but your worksheet has no derivative in it.  Perhaps you meant something else?

## Exact arithmetic...

@Muhammad Usman Floating point approximations can affect eigenvalue calculations severely.  I wouldn't recommend changing exact values to floating point approximations merely to speed up the calculations.  Running your worksheet without the evalfs with M=50 takes less than 8 minutes on my (very old) computer.  About half of that time is spent on calculating the matrix B.  The other half is spent on calculating the determinant.  Solving  the determinant for the eigenvalues takes a fraction of  a second; see the attached worksheet.

8 minutes is not quick, but it's worth the wait for having the assurance that the eigenvalues are computed correctly.

Worksheet: eigenvalues-50.mw

## Why evalf?...

@Muhammad Usman Your MatrixEquation1 is a 20x20 matrix of integers and involves an unknown parameter lambda.  You want to determine lambda so that the determinant of the matrix is zero.  So all you need to do is:

```Z := Determinant(evalm(MatrixEquation1)):
fsolve(Z);```

That produces the 8 roots that you have already obtained.

For some reason that I don't understand, you pollute your equations with several unnecessary evalfs. That introduces numerical truncation errors that lead to the issue that you have observed.  Just don't do evalf and you will be okay.  You don't even need to change Digits.  The default Digits:=10 works just fine.

## Worksheet with corrections...

@apronya Calculate the real and imaginary parts of the solution separately.

 > restart;
 >

 >

 >

 >

 >

 > y0 := R -> cos(w*R);  # find the real part of the solution

 >

 >

 > ds3(parameters=[3]);

 > plots:-odeplot(ds3);

 >