Preben Alsholm

13728 Reputation

22 Badges

20 years, 247 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Christopher2222 Removed spam titles used to disappear quickly, which was helpful: you could see what more you had to remove.
With the new MaplePrimes they don't disappear quickly. This has been mentioned before, but it is really a nuisance. If you click on the title of one of those removed items you are told that the resource if no longer available, and that may very well be because you just removed it yourself. Who wants to remember those silly http titles or what have you? Who is going to use time to make a note of what you have already removed?

@amrramadaneg Let me illustrate by using a system consisting of only two pdes.
This system has the same problem as yours, though. I remedy the situation with a piecewise as in your case.

 

restart; 
sys:={diff(R(x, t), t)=-diff(U(x, t), x)+R(x, t), diff(U(x, t), t)=-(diff(U(x, t), x))};
pdsolve(sys,{R(0,t)=0,U(0,t)=0,R(x,0)=0,U(x,0)=sin(Pi*x)},numeric,time=t,range=0..1); #ERROR as in your case
sys1:={diff(R(x, t), t)=-diff(U(x, t), x)+R(x, t)+piecewise(x>1,1,0)*diff(R(x,t),x), diff(U(x, t), t)=-(diff(U(x, t), x))}; #The new sys.
simplify(sys1) assuming x<=1; #It simplifies to sys, but don't use this below.
res1:=pdsolve(sys1,{R(0,t)=0,U(0,t)=0,R(x,0)=0,U(x,0)=sin(Pi*x)},numeric,time=t,range=0..1);
res1:-plot3d(R,t=2);
res1:-plot3d(U,t=2);
res1:-animate(U,t=1);
res1:-animate(R,t=2);
res1:-value([R,U],t=.5)(.567);
res1:-value([R,U])(0.567,0.5);
res1:-value(R)(0.567,0.5);

See more in the help page:
?pdsolve,numeric

 

@amrramadaneg I tried changing the troublesome equation into

eta*(diff(R(x, t), t))+diff(U(x, t), x)-R(x, t)+piecewise(x>1,1,0)*diff(R(x,t),x) = 0;

where the point is that for x<=1 the last term is just zero. The number 1 was chosen since I used the range 0..1. 
This fools the check in  `pdsolve/numeric/par_hyp` I mentioned in my answer, and a "result" is returned. Whether that "result" is total nonsense I shall let you decide.

It appears that you are using code that I wrote in response to this question of yours:
http://mapleprimes.com/questions/219344-How-To-Solve-Odeproblem#answer232922

I find it appalling that you don't give a link to that yourself. If you had, people would also have a somewhat nicer version to look at.

@shakuntala You ask Tom Leslie about what appears to me my code (from what I can see in this rather confusing forrest of letters without any line breaks).

I find that odd.

If I copy your code and execute, I get a syntax error.
Is this what you intended:
solve(subs(m=ha,f(m))*subs(m=subs(c(t)=a(t),ha), f(m)) = subs(m=ha+subs(c(t)=a(t),ha), f(m)), f);
i.e. the equation to be solved is:

Obviously, if you could solve f(a+c)=f(a)*f(c) for f then you have a solution. You cannot expect solve to do that.
We all know that f=exp is a solution, though.

@chomchom The code follows:
 

restart;
schro := {diff(psi(x), x, x)-(alpha*x^4+x^2-energy)*psi(x) = 0};
ic := {psi(3) = 0, (D(psi))(3) = 1};
Ic := [{psi(3) = 0, (D(psi))(3) = -1}, %$2];
E := [1.06538, 3.306, 5.74796];
schro1 := [seq(subs(energy = e, alpha = .1, schro), e = E)];
##
soln1 := [seq(dsolve(schro1[i] union Ic[i], {psi(x)}, numeric,output=listprocedure), i = 1 .. nops(E))]; #output=listprocedure is used instead of default.
##
for i from 1 to nops(E) do N[i]:=evalf(Int(subs(soln1[i],psi(x))(x)^2,x=-3..3,epsilon=1e-8))^(1/2) end do;
##
with(plots):
display(seq(odeplot(soln1[i], [x, psi(x)/N[i]], -3 .. 3, color = [red, blue, green][i]), i = 1 .. nops(E)));
## And extending to -7..7 (without the normalization) as I commented on later:
##
display(seq(odeplot(soln1[i], [x, psi(x)], -7 .. 7, color = [red, blue, green][i]), i = 1 .. nops(E)));

My guess is that in your quantum mechanical model the Hilbert space is L2(-3..3) and not L2(-infinity..infinity).
## As it turns out dsolve can find the solutions without resorting to numerics.
Here the first one as an example.
sol:=dsolve(schro1[1] union Ic[1]);
nm:=evalf(Int(rhs(sol)^2,x=-3..3))^(1/2); #Using numerical integration for the norm.
plot(rhs(sol)/nm,x=-3..3);
## I'm not necessarily implying that you should avoid the numerical approach since the expression in rhs(sol) is rather complicated.

This extrapolation problem to which I have given an answer below, reminds me that I never got any reaction from you to my answer to your question here:

http://www.mapleprimes.com/questions/219344-How-To-Solve-Odeproblem

I would like to hear your reaction to the answer I gave there.

@9009134 In my recollection pdf-files should be no problem.

Here is a test file. It has only one line.
test.pdf

Notice that when you try to upload, you are told that valid extensions are:
doc, xls, txt, zip, pdf, jpg, jpeg, gif, png, mw, mws, mla, ppt, html, htm, msim, docx, xlsx, pptx, qu

@Mac Dude Yes, that arctan(0,0) = 0 follows from the definition of argument(0) as 0, which is written in the help page for argument (here z:=x+I*y):

"For the case where y = 0 , argument(z) = 0 if  0 <= x."

 

@Carl Love Yes, I just noticed. My first deletion after the new interface appeared (you guys must have been busy deleting).
When you don't see that the deleted item is gone and if there are several other items basically looking the same this could be annoying.

@9009134 I try to keep all my dealings with MaplePrimes matters entirely in this forum. Furthermore, if the reference material you have take up that much space I'm afraid I wouldn't take the time to read it or even just browse it.

@loramina123 You give us h(0)=0 D(h)(0)=tehta (theta, I suppose), D(D(h))(infinity)=1/R and say that infinity might be taken as 2 or 3 micron. By a micron I suppose that (in this context with no units) you mean 1e-6.
But then you also say that you have h(2e-6)=200e-9.  But if infinity can be taken as e.g. 2e-6 then you have a conflict of enormous proportions since R = 1e-7, so 1/R=1e7.
Finally you mention D(D(H)(-infinity)=0 and that -infinity can be equal to 2e-6. I suppose that you meant -2e-6?

In these matters one has to be more careful or else people give up.

@chomchom Are you sure that you want to integrate from -infinity to infinity?

If you are, then try e.g.

display(seq(odeplot(soln1[i], [x, psi(x)], -7 .. 7, color = [red, blue, green][i]), i = 1 .. nops(E)))

If you think that is fine then replace -7..7 by -19 .. 19.

@loramina123 So what are the boundary conditions? You have presented more than one set of conditions. The ode appears to be the same though.
The original 3-point boundary condition problem may be doable, but you do have to use 4 conditions for a 4th order ode. So two conditions could be given at 0, e.g. values of h(0) and D(h)(0). Besides one condition at each of the other two points. By using smoothness at the middle point this can (sometimes) be done.
I did as a rather trivial example the following, where an exact solution can also be found. Thus we can compare the numerical result with the exact.

ode:=diff(y(x),x$4)+4*y(x)=0;
sol:=dsolve({ode,y(0)=0,D(y)(0)=1,y(Pi/2)=0,y(Pi)=0}); #This gives you a simple formula for the solution

By doing what I sketched above I got very good agreement between numerical and exact results.

Before attempting anything similar to your problem I need the actual boundary conditions if you are still interested in the 3-point problem.

First 74 75 76 77 78 79 80 Last Page 76 of 230