Preben Alsholm

13743 Reputation

22 Badges

20 years, 340 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@iman Well, actually I expected you to do the work, not me.
You may want to read the help for D.
?D
and the help for pdsolve/numeric:
?pdsolve,numeric
See the examples.

@golnaz You are right. The inverse of the incomplete beta function evaluated at half the argument is not what I gave.
There is still confusion about the square root, but that confusion has nothing to do with that issue:

restart;
IB(x,a,b)=Int(t^(a-1)*(1-t)^(b-1),t=0..x); #This is the definition of the incomplete beta function
##Now we find it as a function of x,a,b:
B:=int(t^(a-1)*(1-t)^(b-1),t=0..x);
IB:=unapply(B,x,a,b);
##The confusion mentioned is: Did you mean the following one instead?
#IBsqrt:=(t,a,b)->IB(sqrt(t),a,b);
plot(IB(x,1/5,1/2),x=0..1); #The graph of IB as a function of x.
plot([IB(x,1/5,1/2),x,x=0..1]); #Graph of inverse of x->IB(x,1/5,1/2) (call it invIB )
plot([2*IB(z,1/5,1/2),z,z=0..1]); #This is the graph of invIB(f/2) as a function of f.

The argument is that the graph of invIB(f/2) we can get as
[f, invIB(f/2), f=0..2*IB(1)]
This graph is also obtained by using z=invIB(f/2), i.e. f=2*IB(z):
[2*IB(z),z,z=0..1]
##I edited the answer above to correct this error. Notice though that there I still use the square root.


I'm somewhat surprised (to put it mildly) that you don't link to your previous question:

http://www.mapleprimes.com/questions/205318-Error-In-DEplot#comment219317

@rit I used the right hand sides 0, 3*t^2, and 2*t^3. Notice that I was in doubt about what you meant them to be. Clearly with these right hand sides neither system is autonomous, so the concept of sink or source is not applicable.

So give us the system with equality signs, so that the question "autonomous or not?" can be settled. To find equilibria and determine eigenvalues for the jacobian at those points can then be done.

You have three unknowns (a(t), b(t), c(t) ) and  two odes.

I copied your 2 lines and had no problem. The result returned was -36949.67738.

Maybe you should try your own 2 lines again after a restart.

@emmantop You copied my code incorrectly. It should be

dsolve(eval({ode1, ode2} union BCS, E union {k[1]= 0}), numeric, method = bvp[midrich], abserr = 1.10^(-8))

Notice that there we have {k[1]=0} and NOT {k[1]}=0.

You had
`union`(E, {k[1]}) = 0

which is equivalent to

E union {k[1]} = 0

which is (basically) nonsense, but not unsyntactical (in Maple).

Let me add that what I referred to as "weird" is really just the zero solution for f. When D(f)(0)=1 is removed then all f related boundary conditions are homogeneous and f(eta)=0 solves ode1.
Another thing while at it:  Did you really want abserr = 1.10^(-8)? That would be abserr=0.4665073802.
I had success (i.e. in two cases as before mentioned) with abserr=1e-5.

@Axel Vogt Epsilon prints like a capital epsilon, which ought to be fine.

By the way I don't see the solution found by Carl in your reduced system.

Another thing: I cannot make fsolve give a result on Sys[2] without omitting convert/rational. Without that, no problem.

Using solve on the equations without assignment to the constants gives two real solutions altogether when assigning afterwards:

For ppsi = .1561816797:
{A1 = .2178679378, A2 = 0, C0 = 0, C1 = 0, C2 = 0, E0 = 2.923409461, H = 1.169518511, H1 = -3.370336272, K = .2120924177, L0 = -1.923409459, L1 = 1, R = 0.1077611106e-1, W = .4031652600, Y = .6858871836, Epsilon1 = -1.108403687}
and
{A1 = .1151377592, A2 = -0.2771025818e-1, C0 = -0.2106664520e-1, C1 = -0.3108676424e-1, C2 = -0.4587284305e-1, E0 = 2.923409461, H = .9700004405, H1 = -.4209490949, K = 0.8582496732e-1, L0 = -3.093516599, L1 = -.1854153076, R = .655446246, W = .3221564949, Y = .4545705612, Epsilon1 = .2686838282}

For ppsi=1.947717707:
{A1 = 0.1146472975e-1, A2 = 0, C0 = 0, C1 = 0, C2 = 0, E0 = .2344194943, H = 0.9509540701e-1, H1 = 0.6343537741e-1, K = 0.1116080812e-1, L0 = .7655805057, L1 = 1, R = .3632312348, W = .3518971259, Y = 0.4867856263e-1, Epsilon1 = 0.3694572998e-1}
and
{A1 = 0.1429702382e-2, A2 = 0.5359982543e-1, C0 = 0.4074911503e-1, C1 = 0.6013098524e-1, C2 = 0.8873162983e-1, E0 = .2344194943, H = .5898272830, H1 = 1.456583019, K = 0.5218750955e-1, L0 = .1114825845, L1 = .3373447074, R = .655446246, W = .3221564946, Y = .2764103064, Epsilon1 = 0.4295105618e-1}

The above results were computed with Digits=10, but confirmed with Digits=15.

(NOTE: Second value of ppsi).
Since you know the solution you can use that as a starting point for fsolve just to verify the result, always a good idea.
With the method I mentioned in my previous comment (but didn't elaborate on) I found

{A1 = 0.142970229228308e-2, A2 = 0.535998252919882e-1, C0 = 0.407491149283945e-1, C1 = 0.601309851312542e-1, C2 = 0.887316296122945e-1, E0 = .234419494359216, H = .589827281889861, H1 = 1.45658301648039, K = 0.521875094397844e-1, L0 = .111482584275646, L1 = .337344707782075, R = .655446246865986, W = .322156494043605, Y = .276410305860281, Epsilon1 = 0.429510561240001e-1}

which is the one you are talking about.

Using a starting point can help fsolve quite a lot. So any method that will provide you with an approximate solution may help.

I was using your system as a test on some solving method I'm looking into and found one more solution (there may be more):
(NOTE: First value of ppsi).
It has  
sol1:={A2=0,C0=0,C1=0,C2=0,L1=1};
#and
sol2:={A1 = .217867938869731, E0 = 2.92340946074119, H = 1.16951851156081, H1 = -3.37033625373258, K = .212092418616612, L0 = -1.92340946074119, R = 0.107761071012707e-1, W = .403165260852399, Y = .685887185020117, Epsilon1 = -1.10840368196766};

#You can test the soluion with fsolve by doing

eval(eq,{A2=0,C0=0,C1=0,C2=0,L1=1});
neweq:=select(hastype,%,name);
fsolve(neweq,sol2);
#However, the following gave me the solution found by Carl (I'm also using Digits=15):

fsolve(eq, sol1 union sol2);

@9009134 I'm looking at your latest worksheet called frekans.mw. I don't know what you call the original.

1. I don't see any code like that in that worksheet.
2. The try ... end try construction is used when you want to try something, but expect that some error might occur and don't want that to stop the process from running. Here you want to try different extra bcs.
If dsolve doesn't work with that extra bcs (b) and results in an error then you want that error caught (the catch clause). In case of error the last error is printed, but the loop is not interrupted, you go on to the next b in extra_bcs.
3. Indeed! Try changing it.
4. b is a loop variable just as in this simple construction:
   for b in [7,9,13] do b^2 end do;
4 II. Only the b's that didn't result in errors produce results, so you can pick any of those.
5. Since RES was defined as res[b] for one of the successful extra bcs b, RES(0) will give you the value of the solution for theta=0. That was just intended for getting the value of omega, which supposedly is constant. For that worksheet I get 9.03506039191669.
You might as well have done RES(1);

@9009134 dsys4 uses theta. You use q later. Change each q to theta (or vice versa).

@Kitonum SCR stands for "Software Change Request" mostly thought of as a euphemism for bug report.
You can find a way to do it in the menu line on top on this site under "More".

You start with a substitution into Ca[2] and receive an answer involving sin. Where did that come from?

If we are to help you we must see the code from the start. You can forget about the output: remove it before copying.

@Al86 Just do:

eval~(Ymono,t=~T);

First 107 108 109 110 111 112 113 Last Page 109 of 231