Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@mobiusinfi You can only have 2 boundary conditions. So think again.
I tried:

restart;
k :=21.37595739:
c :=1.760652758:

v:=unapply(u(x,t)+k*u(x,t)^(2),x,t);

PDE := diff(u(x, t), t)-(diff(v(x, t), x, x)) = 0;
#Leaving one of your boundary conditions out:
IBC1 := {u(x, 0) = 0, D[1](v)(0, t) = 0, D[1](v)(1, t) = c};
S1 := pdsolve(PDE, IBC1, numeric, time = t);
S1:-animate(t=5);
#Leaving another out instead produces the error known from earlier:
IBC2 := {u(1,t)=1,u(x, 0) = 0, D[1](v)(0, t) = 0};
S2 := pdsolve(PDE, IBC2, numeric, time = t);
S2:-animate(t=5);


@tomleslie I doubt that the LambertW solution is going to solve the problem. I tried the following:


restart;
PDE := diff(u(x, t), t)-diff(u(x, t)+k*u(x, t)^2, x, x) = 0;
sol:=pdsolve(PDE);
#This solution appears to be of the form
sol2:=u(x,t)=A*LambertW(-exp(a*x+b*t+c))+B;
# where A, a,b,c, B are 5 arbitrary constants (not necessarily real!)
#However, their number is reduced to 3:
pdetest(sol2,PDE);
collect(numer(%),LambertW);
coeffs(subs(LambertW(-exp(a*x+b*t+c))=L,%),L);
solve({%},{B,b});
expand(%);
sol3:=eval(sol2,%);
#Clearly sol3 is not going to solve the initial value problem:
eval(sol3,t=0);

What we are getting is more likely just a building block as when pdsolve gives us in case k=0:
pdsolve(eval(PDE,k=0),build);





I looked at your pde in the following form, where your k = 21.37595739 and your c = 1.760652758.
Even using k=0 gave problems!

restart;
PDE := diff(u(x, t), t)-diff(u(x, t)+k*u(x, t)^2, x, x) = 0;
IBC := {u(1, t) = 1, u(x, 0) = 0, (D[1](u))(1, t) = c};

res:=pdsolve(eval(PDE,k=0.00001),eval(IBC,c=1),numeric,time=t,range=0..1);
res:-animate(t=.3);

What kind of behavior are you expecting? Your problem obviously comes from some physical situation. What is that?

Example:

You need the independent variable x in the ode

diff(y(x),x) = y(x)

You cannot just write

diff(y(x),x) = y

That is very likely your problem. There may be more problems.

@patient But your sol is not a solution of sys in your latest worksheet numeric_solution2.mw:

odetest(sol,sys);
simplify(%) assuming z::real,t<0;

You see that the second member is not identically zero.


Now try with coefficients A and B:
sol1 := {u(z) = (1+(1/4)*z)*exp((1/8)*z)*A, v(z) = exp((1/4)*z)*B};
odetest(sol1,sys);
simplify(%) assuming real;
collect(%,z);

Notice that we need A=sqrt(B)  (unless B=0).

So the question is: What do you want B to be?


@patient 
I wrote at the end of my reply Exact solution:
"The initial conditions which this solution satisfies at x00=sqrt(x0*t) are easily computed and there is a discrepancy between those and the ones I named ICS in your numerical approach."

Taking a look at that in the present context:

ICS; #Initial conditions you use for U and V
SOL; #The exact solution for U and V
#We find the initial conditions actually satisfied at x=x0=sqrt(t*z0) by the exact solution:
ICS0 := eval(SOL, x = sqrt(t*z0));
diff(SOL[2], x);
convert(%, D);
eval(%, x = sqrt(t*z0));
ICS1 := simplify(%) assuming t<0;
ICSE :=ICS0 union {ICS1};


You see that there is quite a discrepancy between ICS and ICSE.

@patient x0 (my x00) is to be treated as a parameter, so must be a name. Thus you should leave out the assignment to x0.
Leaving that out everything works.

@patient I had a look at your worksheet called exact_solution.mw.

I suggest shortening what you are doing there to:

restart;
u1 := -(-1-(1/4)*x^2/t)/(-t)^(3/2);
u2 := exp((1/8)*x^2/t);
F:=unapply(piecewise(x^2/t>-4,u1*u2,0),x,t);
plot([seq(F(x,-0.01*i),i=1..10)],x=-1..1,color=black);
##OK the last curve is not red, but in this animation it is:
plots:-animate(plot,[F(x,t),x=-1..1,color=`if`(t<-0.01,black,red)],t=-.1..-0.01,frames=10,trace=9);

So are you claiming that U(x) = F(x,t) solves the system SYS obtained above with some appropriate S(x)?
And what is that S(x) supposed to be?
Why is x^2/t >-4 appearing as a break?

########################
An exact solution can be found from earlier remarks. You confused me by renaming v(z) to s(z).
sol:= {u(z)=(1+z/4)*exp(z/8), s(z)=exp(z/4)};
odetest(sol,sys);
simplify(%) assuming z::real; #OK
##So a solution to SYS is
SOL:=subs(s(z)=S(x),u(z)=U(x),z=x^2/t,sol);
odetest(SOL,SYS);
simplify(%) assuming real;

The initial conditions which this solution satisfies at x00=sqrt(x0*t) are easily computed and there is a discrepancy between those and the ones I named ICS in your numerical approach.



@Kitonum This is somewhat shorter:
simplify(temp) assuming positive;
radnormal(%);


@Christopher2222 Sorry about not catching the obvious mistake right away. Assuming that up is the same as up on your figure in the worksheet, then the altitude is given by h=R*sin(theta(t)) when theta is measured from the positive x-axis as you indicate in your figure.

The potential energy is then m*g*h + C, where the constant C can be chosen arbitrarily.

Had the angle theta been measured from the positive y-axis then the altitude would be given by h=R*cos(theta(t)) and the potential energy would again be m*g*h + C. Keeping the initial conditions then the events trigger would be theta(t)=Pi/2 (and theta(t)=-Pi/2, but that won't ever happen for your initial condition).

@patient Clearly the initial conditions for S and U must be different from the initial conditions for s and u. The change of variable from z to x obviously has an effect.
I won't change your (to my mind) unfortunate name x0 for the initial z-value (z0 appears to me to be more natural).
I shall call the x-value corresponding to z=x0 for x00.
That value is found from x^2=t*z to be x00 = sqrt(t*x0) if we choose x >0, which I shall do.

The derivative of S at x = x00 is found by using the chain rule to be
S'(x) = s'(z)*2*x/t
thus
S'(x00)= s'(x0)*2*x00/t.

Now we notice that the parameter t appears in the initial conditions as S(sqrt(t*x0)) = ...
It turns out that dsolve doesn't accept that, whereas the form S(x00) = ... with x00 as a parameter is accepted. Whether that should be considered a bug or not is irrelevant right now: we can work around it by using x00 as parameter instead of t.

Thus we can do:
SYS2:=subs(t=x00^2/x0,SYS);
ICS:={S(x00)=exp(x0^2-16),U(x00)=9*exp(x0^2-10),D(S)(x00)=-4*exp(-12)*2*x0/x00};
res:=dsolve(SYS2 union ICS, numeric, method = rkf45_dae, parameters = [x00]);
res(parameters=[sqrt(x0*(-.1))]); # This means that t = -0.1
plots:-odeplot(res,[[x,100*S(x)],[x,U(x)]],0..10);

Notice that I have multiplied S by 100 so it doesn't get drowned by U.

The procedure p can be modified as follows:
p := proc (t)
   res(parameters = [sqrt(t*x0)]);
   plots:-odeplot(res, [x, U(x)], 0 .. 4)
end proc;

plots:-animate(p, [t], t = -.1 .. -0.1e-1,view=0..0.05);


@Christopher2222 The zero for potential energy can be set arbitrarily. It is of no consequence for the equations of motion. In fact just look at the case in hand. L = T-U. You only use two derivatives of L. Both would eliminate an additive constant to U. The main thing is that since the gravitational pull is downwards, U should increase with the altitude.

@Christopher2222 Obviously the potential energy should increase with the altitude. Without a change of sign it doesn't. Whether U is positive or negative is irrelevant.

As Rouben points out, there is nothing wrong.

I can produce the same error by moving the local declaration:

game := proc()
  #local player1, player2, roll;
  roll := rand( 1..6 );
  local player1, player2, roll;
  player1 := roll():
  player2 := roll(2):
  if player1>player2 then "A wins"
  elif player1=player2 then "Tie"
  else "B wins"
  end if;
end proc:

I think you ought to reveal some more of what you are doing. I sure don't understand it.

First 120 121 122 123 124 125 126 Last Page 122 of 231