Preben Alsholm

13728 Reputation

22 Badges

20 years, 244 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

We notice in the numerical solution that R doesn't vary much, between 0.14 and 0.1428.
Thus it might be illuminating to look at the special case where R(theta) in the ode for mu is replaced by a constant R0.
It turns out that this ode with the initial conditions mu(0) = mu0, mu'(0) = 0 has a solution satisfying mu'(Pi) = 0 provided only that R0 and mu0 satisfy the inequality 
-mu0-(1/2)*ln(R0)+(1/2)*ln(2) < 0.
 

restart;
ode0:=R0 = 2*exp(-2*mu(theta))*(1-(diff(mu(theta), theta, theta))-cot(theta)*(diff(mu(theta), theta)));
res:=dsolve(ode0);
resA:=allvalues(res);
## The second result equals the first if _C1 is replaced by -_C1:
resA[1]-eval(resA[2],_C1=-_C1);
## Thus we need only consider the first result resA[1].
###### Graphical experimentation:
eval(rhs(resA[1]),{R0=0.14,_C2= -1,_C1=1});
plot(eval(rhs(resA[1]),{R0=0.14,_C2=-.1,_C1=1}),theta=0..Pi);
plots:-animate(plot,[eval(rhs(resA[1]),{R0=0.14,_C2=-2}),theta=0..Pi],_C1=.9..1.1);
###### Now look at mu'(theta):
diff(rhs(resA[1]),theta);
eval(%,{_C1=1}); #Seems to be a good idea
limit(%,theta=0,right) assuming _C2<0,R0>0;
limit(%%,theta=Pi,left) assuming _C2<0,R0>0;
## So the following expression satisfies mu'(0) = mu'(Pi) = 0 as long as _C2<0 and R0>0:
sol1:=eval(rhs(resA[1]),_C1=1);
## We want mu(0) = mu0 (at least as a limit from the right):
limit(sol1,theta=0,right) assuming _C2<0,R0>0;
c2:=solve(%=mu0,{_C2}); # using that mu(0)=mu0
sol:=eval(sol1,c2); ## Solution for general R0 and mu0
For _C2 to be negative we need R0 > 1/5000:
ineq:=rhs(op(c2))<0;
## The parameters in your case:
params:={R0=0.1428,mu0=ln(100)}; # The "final" value from numerics is taken for R0.
evalf(eval(ineq,params)); #OK satisfied
eval(sol,params); #See the actual solution
plot(eval(sol,params),theta=0..Pi,labels=[theta,mu]); 
plot(eval(diff(sol,theta),params),theta=0..Pitheta=0..Pi,thickness=3,size=[800,default],labels=[theta,diff(mu(theta),theta)]) 

 

I get no errors in that whole worksheet.
interface(version);
  `Standard Worksheet Interface, Maple 2017.3, Windows 10, September 27 2017 Build ID 1265877`
kernelopts(version);
  `Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877`

@Thomas Richard It appears to me that the `*` used and vartheta instead of theta is not a problem affecting the solution at all.

So what in your opinion should C, I, and R be?

There is no attachment!

I don't see any attachment.

@vv Just a comment on infolevel: With  infolevel[Basis]:=3:  (or 2) we do get information, but it may not be of interest to codell.
Since the procedure remembers the information will only come the first time it is run with the given input.

 

I have never heard of acximetric tolerance. Is it spelled correctly?
You could be a little more specific and provide an actual example.
When you say " For example +10/-5" I still don't get what you are talking about.

I read at the top of the pdf file containing the problems as formulated by your teacher:
"Studentene får hver sin oppgave, som skal løses selvstendig."'

Does that not exclude getting help from MaplePrimes? Ask your teacher.

@tomleslie I have noticed before that these problems come up too often for students in Denmark.
I cannot read the uploaded file either.

Personally, I haven't had these problems in years although my computer is made for use with Danish characters and bought in Denmark.
Years ago when I occasionally had problems, I started a habit of never (well, almost never) saving a worksheep with output.
I routinely remove all output from the worksheet before saving it by going to
Edit/Remove Output/From Worksheet.
Another thing is that I never use 2D math input or fancy text. My worksheets are rather plain, but I do use the Standard Interface (with Maple Input and Worksheet mode).
##
acer has been able to rescue some people in the past.

@Markiyan Hirnyk Where do you find HINT=sum documented?
HINT=xxx works the same:
 

restart;
PDE := (y-u(x, y))*diff(u(x, y), x)+(u(x, y)-x)*diff(u(x, y), y) = x-y;
infolevel[pdsolve] := 3: 
ans := pdsolve(PDE, HINT = xxx);

 

@oceanxinlie To answer your question: 
Why do you set "approxsoln=[y(x)=tanh(x)]" ?
let me first say that dsolve/numeric/bvp has to construct a starting solution itself (an approximate solution) if none is given by the user.
To do that it uses the boundary conditions only, i.e. y(0) = 0 and y'(0) = 0. It uses a polynomial of lowest possible degree satisfying those conditions. In this case that is the trivial y(x) = 0.
Since we might have some idea about the form of the solution it is often a good idea to give an approximate solution. It need not be very good. In this case the one I used happens to be rather bad since tanh(x) > 0 for all x > 0.
But it happened to work!
Among nonzero polynomial approximate solutions this one is an example and it works fine:
y(x) = (1/2945)*x^2-2*x.

@oceanxinlie I was a little slow in seeing the obvious: The solution to your ode with the boundary conditions as given can never attained a positive value.
Here is the proof:
First of all notice that if the ode is written as
y''(x) - a*y(x)^0.337 = b*x, a > 0, b > 0,
then for the term y(x)^0.337 to make sense (be real) we must demand y(x) >= 0 for all x.
But from the ode it follows that
y''(x) = a*y(x)^0.337 + b*x
thus y''(x) >= b*x, so by integrating from 0 to x we get y'(x) - y'(0) >= (b/2)*x^2.
But y'(0) cannot be negative since then y(x) would be negative for some x near 0.
Thus it follows that y'(x) >= y'(0) + (b/2)*x^2  >= (b/2)*x^2 > 0 for x > 0.
So the solution cannot satisfy y'(2945) = 0.
## Suppose then that the ode is rewritten as y''(x) - a*abs(y(x))^0.337 = b*x.
Then there it isn't different from the first one if y(x) >= 0 all x. Thus the same conclusion.
Now suppose that there is an x1 with y(x1) < 0. If the solution is to assume a positive value for some x2 > x1 then we may assume that also y'(x1) > 0 since if no such x1 < x2 existed y(x2) couldn't be positive.
Thus we have y''(x) = a*abs(y(x))^0.337 + b*x >= b*x.
So by integrating from x1 to x > x1 we get
y'(x) - y'(x1) >= (b/2)*(x-x1)^2
and therefore
y'(x) >= y'(x1) + (b/2)*(x-x1)^2  >= (b/2)*(x-x1)^2 > 0 for x > x1,
which makes y'(2945) = 0 impossible.

Using another revised version still gives a negative result and one that is about 1.75*10^5 times smaller in absolute value than the one where abs is used:

restart;
Digits:=15:
u:=piecewise(y(x)>=0,y(x)^0.337,0);
ode2:=diff(y(x),x,x)-0.00003019*u=9.542*10^(-13)*x;
bcs:=y(0)=0,D(y)(2945)=0;
res2:=dsolve({ode2,bcs},numeric,method=bvp[midrich],approxsoln=[y(x)=tanh(x)]);
plots:-odeplot(res2,[x,y(x)]);

## NOTE: This 'desperate' attempt to get a positive solution by rewriting the ode ends up just giving the same as the one where u = 0.

@Christopher2222 Yes, in my experience you can set libname only once in a session and expect it to work.

First 56 57 58 59 60 61 62 Last Page 58 of 230