Preben Alsholm

13728 Reputation

22 Badges

20 years, 252 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@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);

@rit

I don't know what you mean by "I choose one or zero". One or zero where?

The RootOf result from solve will give 3 different sets of solutions for (a,b,c). Choose the one that gives you real solutions in a neignorhood of t=1 (if real solutions are what you want).
The other two solutions will be imaginary.
Forget about symbolic solutions.
Example:

eval[recurse](convert({ode},D),{t=1,D(b)(1)=0}); #Example
solve(%,{a(1),b(1),c(1)});
## So it looks like we can get away with:
ics1:=a(1)=1,b(1)=0,c(1)=1; #Example
res:=solve({ode},diff({a,b,c}(t),t));
Ares:=allvalues(res):
nops([Ares]);
Ares[1];
sol1:=dsolve(Ares[1] union {a(1)=1,b(1)=1,c(1)=1},numeric);
plots:-odeplot(sol1,[[t,a(t)],[t,b(t)],[t,c(t)]],0.5..1.62);
plots:-odeplot(sol1,[t,a(t)],0.5..1.62,thickness=2); p1:=%:
sol2:=dsolve(Ares[2] union {a(1)=1,b(1)=1,c(1)=1},numeric);
sol3:=dsolve(Ares[3] union {a(1)=1,b(1)=1,c(1)=1},numeric);
sol2(1.2345);
sol3(1.2345);
plots:-odeplot(sol2,[[t,Re(a(t))],[t,Im(a(t))]],-1..5,color=[green,red]);p2:=%:
plots:-odeplot(sol3,[[t,Re(a(t))],[t,Im(a(t))]],-1..5,color=[green,red]); p3:=%:
plots:-display(p1,p2,p3);




@taro The level of evaluation inside a procedure has nothing to do with the fact that your procedure g:=x->f1 doesn't work.
The reason is that before any evaluation at all takes place (no matter the level) when you do g(2); the variable x used in x->....  is replaced by 2 (simple substitution). Now no x is seen in the body. It will be seen after evaluation surely, but that will take place after the substitution and that is too late, so that x doesn't get the value 2.

I tried using the first part of your code till you have defined all the equations and the variables, eqs and vars, respectively.

Then I used fsolve:

res_phi:=fsolve(eqs,{vars[]}=~1); #Giving start values 1 to all variables

After that I defined the matrix M by

M:=Matrix(subs(res_phi,[seq([seq(phi[p,q],p=0..Numy+1)],q=0..Numz+1)]));

This matrix I compared to your matrix sol3. The difference matrix M-sol3 had elements of absolute value less than 5e-11:

max(abs~(M-sol3));
     4.98876495669264841*10^(-11)

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