Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

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)

@Carl Love As for second order in time we can take the wave equation, which of course is not elliptic, but you said "any PDE that isn't first-order in time".

We can compare with the exact solution. Incidentally the exact solution found by pdsolve is wrong!

restart;
pde:=diff(u(x,t),t,t)=diff(u(x,t),x,x);
ibcs:=u(x,0)=0,D[2](u)(x,0)=sin(x),u(0,t)=0,u(Pi,t)=0;
sol:=pdsolve({pde,ibcs}); #WRONG! But we can see that u(x,t) = sin(x)*sin(t) is a solution.
pdetest(u(x,t)=sin(x)*sin(t),[pde,ibcs]); #Check
eval((rhs-lhs)~({pde,ibcs}),u=((x,t)->sin(x)*sin(t))); #Double check
res:=pdsolve(pde,{ibcs},numeric); #Numerical solution
res:-animate(t=4*Pi,frames=100); p1:=%:
plots:-animate(plot,[sin(x)*sin(t),x=0..Pi],t=0..4*Pi,frames=100); p2:=%:
plots:-display(p1,p2); #Doesn't look bad at all!


@Carl Love I try, but often fail, to choose my wording carefully.
I wrote
"Thus the IBC given by
IBC2:={u(0, y) = 1, D[1](u)(1, y) = 0, u(x, 0) = 1, D[2](u)(x, 0) = 0};
would not be rejected by pdsolve/numeric." (emphasis added).
That was what I meant.

##############
It may be fair to say that the time-based solver is not doing a great job on elliptic pdes, but it isn't failing entirely.

Take the laplace equation:

restart;
pde := diff(u(x, y), x, x)+diff(u(x, y), y, y) = 0;
U:=evalc(Re(exp(x+I*y))); #Harmonic, so solves pde
#u(x,y) = U satisfies
ibcs:={u(0, y) = cos(y), D[1](u)(1, y) = exp(1)*cos(y), u(x, 0) = exp(x), D[2](u)(x, 0) = 0};
pds:= pdsolve(pde,ibcs, numeric,timestep=0.001);
pds:-plot3d(y=0..0.3); p1:=%:  #Yes, don't go too far with y. Already a ripple starting.
plot3d(U,x=0..1,y=0..0.3); p2:=%:
plots:-display(p1,p2);
pds:-value()(.2345,0.07913);
eval(U,{x=.2345,y=0.07913});


Let me add that the error message may be misleading. It is not so much a question whether the pde itself is elliptic or not. You have to include the boundary conditions, IBC in your case. They are given on the square [0,1]x[0,1]. Thus the problem is not handled by the algorithm used by pdsolve/numeric.
This algorithm is time based, which basically means that it starts from time t=t0 with enough information to march forward step by step. In your case time could be taken as y. Since the pde is second order in time (y) in order to get going the values of u(x,y0) and D[2](u)(x,y0) should be known (for some y0,here 0, and all x in some interval, here [0,1]).

Thus the IBC given by
IBC2:={u(0, y) = 1, D[1](u)(1, y) = 0, u(x, 0) = 1, D[2](u)(x, 0) = 0};
would not be rejected by pdsolve/numeric. You can give the optional argument 'time'=y, but it is not necessary since it can be deduced from IBC2.

First you have to correct a few syntax errors: missing arguments in 'u' in two places, a space between D[i](u) and what follows in the IBC (disastrous in 2D-math input: avoid it like the plague!).

After that you got:

PDE := diff(u(x, y), x, x)+diff(u(x, y), y, y) = u(x, y)^(1/2)+diff(u(x, y), x)^2/u(x, y)^(3/2);
IBC := {u(0, y) = 1, u(x, 0) = 1, D[1](u)(1, y) = 0, D[2](u)(x, 1) = 0};
pds := pdsolve(PDE, IBC, type = numeric);

The result is the error message:

Error, (in pdsolve/numeric) unable to handle elliptic PDEs



@necron Now with names (and more) changed I again suggest that before attempting to plot you should try to see if you can get values out of quantities like

evalf(g2(2.345)); No chance if you cannot do these:
evalf(eval(c2,{u=1,r=2})); #Comes up with a number
evalf(eval(c2,{u=1,r=20})); #This doesn't
evalf(eval(c3,{u=1,r=2})); #This doesn't either


@Prashanth I think this little session shows what happens. I'm using the same example as I did above.

restart; #Important
p := proc()
    with( ListTools );
    Reverse( [ 1, 2, 3 ] )
end proc:
eval(Reverse); #with(ListTools) has not been executed, but Reverse inside is now the global :-Reverse
p(); #This will cause execution of with(ListTools)

eval(Reverse); #This is ListTools:-Reverse
Reverse([1,2,3]); #It works
p(); #Reverse inside p remains :-Reverse
op(0,p()); # This is :-Reverse
%-Reverse; # Not zero
###Now executing without a restart the same procedure definition once more
p := proc()
    with( ListTools );
    Reverse( [ 1, 2, 3 ] )
end proc:
eval(Reverse); #with(ListTools) was already executed above so the Reverse seen is ListTools:-Reverse
p();
###If you try the whole thing from the start and use another name (q) for the procedure the second time, then you will have both p and q. p still doesn't work, but q does.



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