Preben Alsholm

13728 Reputation

22 Badges

20 years, 242 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

Your domain is B = { (s,t) | s in [0,L] and t >=0}.

The gradient of w(s,t) is the vector < diff(w(s,t),s), diff(w(s,t),t) >.
A unit normal to the part of the boundary given by (s,t) = (L,t), t>=0 is v = <1,0>.
The normal derivative is the dot product of the gradient and v, thus equal to diff(w(s,t),s).

However, your IBC-conditions contain diff(w(L,t),t) on the boundary s=L.

You can use Flatten from ListTools and then nops. If you mean that no number should be counted twice, you can use MakeUnique.

S:=[seq([seq(i, j = i .. 5)], i = 1 .. 5)];
FS:=ListTools:-Flatten(S);
nops(FS);
ListTools:-MakeUnique(FS);
nops(%);
#Obviously, if you start with a flat version you don't need Flatten:

S2:=[seq(seq(i, j = i .. 5), i = 1 .. 5)];
nops(S2);
nops(ListTools:-MakeUnique(S2));

@newlam We haven't seen your corrected version, but if you corrected the occurrence of wx(L,0)=0 to wx(L,t)=0, then it follows that D[2](wx)(L,t)=0.

There appears an H without value in IBCX and IBCY. Also there is a ws in IBCZ. Presumably these are typos.

First try the following with Digits = 10 (default), then try setting Digits to e.g. 20.

After that try an exact approach i.e. replace floating point numbers by rationals.

restart;
Digits:=20:
eqs := {2.32 = a+b, pa = .4310344828*a, pb = .4310344828*b, pa+pb = 1};
solve(eqs);
T:=LinearAlgebra:-GenerateMatrix(eqs,[a,b,pa,pb],augmented);
LinearAlgebra:-GaussianElimination(T);
LinearAlgebra:-LinearSolve(T);

restart;
eqs := convert({2.32 = a+b, pa = .4310344828*a, pb = .4310344828*b, pa+pb = 1},rational,exact);
solve(eqs);
T:=LinearAlgebra:-GenerateMatrix(eqs,[a,b,pa,pb],augmented);
LinearAlgebra:-GaussianElimination(T);
LinearAlgebra:-LinearSolve(T);

Added. The problem can be reformulated

a + b = c
s*a+s*b = 1
where s = .4310344828 and c = 2.32

This system has no solution unless 1-s*c = 0, in which case it has infinitely many.
Now s*c =1.000000000 with Digits set at 10, but 1.0000000001 with Digits set at 11, so the question is one of round-off.

In the assignment to PDEZ there appears the name wz. It ought to have been wz(s,t).

Try

indets(PDE,name);
                           {s, t, wz}
indets~([PDEX,PDEY,PDEZ],name);
                  [{s, t}, {s, t}, {s, t, wz}]

and try to run the worksheet without 'declare(....)' to find the error.

Have you tried the frames option?

The help page for odeplot:

"If the frames=n option is given, then odeplot produces an animation of the solution over the independent variable range for the plot, from lowest to highest value. The final frame of the animation is the same as the plot created without the frames option. Note: The refine option cannot be used with animations, as odeplot must choose the points for the animated plot."

You can use dchange from the PDEtools package.

The change of variables used here is  x=-2+alpha,y(x)=-1+v(alpha)

ODE := diff(y(x), x) = (y(x)+1)/(x+2)-exp((y(x)+1)/(x+2));

PDEtools:-dchange({x=-2+alpha,y(x)=-1+v(alpha)},ODE,[alpha,v]);

#You can do the whole reduction to a separable equation like this:

PDEtools:-dchange({x=-2+alpha,y(x)=-1+alpha*u(alpha)},ODE,[alpha,u]);
isolate(%,diff(u(alpha),alpha));

It seems that I ought to have seen this problem handled before, but maybe not.

IntegrationTools has an Expand that does what you expect

IntegrationTools:-Expand(Int(f(x)+g(x),x=a..b));
            
However, I don't recall seeing a similar thing for Sum (or sum).

It is easily written though. Something like this.

ES:=proc(S::specfunc(anything,{sum,Sum})) if type(op(1,S),`+`) then map(op(0,S),op(1,S),op(2,S)) end if end proc;

Then the following commands result in zero

S:=sum(sum((b[i]-b[j]),i=1..n),j=1..n);
ES(S);
simplify(%);




I don't think you will get any symbolic answer.

Replacing 'int' by 'Int' even

s1:=Int(.......)
eval(s1,{k=1,sigma=1,Q=0,theta=0});

#is not found (in Maple 15)

value(%);

#But of course you can get a numerically computed result

evalf(%);
                         0.01405865281



You are using curly braces { } where you should use square brackets [ ] so that your intended order is respected.

Try this line instead where {A,B} is replaced by [A,B].

AAA := LinearAlgebra[Eigenvectors](op(subs(s1 = 0.19e-3, [A, B]))):

It is not surprising that you get two different results when using [A,B] and [B,A].

Notice that if f is a solution to the boundary value problem, so is -f. So if there is a solution, then it is not unique. It seems that there are in fact at least 4 solutions.

You must solve the problem numerically. I didn't have any luck using dsolve/numeric directly on your problem, so I turned your problem into an initial value problem with the parameter p being D(f)(-lambda).

We must determine p s.t. f(1) = 0.

(I used the values given for beta and lambda).

init:=dsolve({s,f(-lambda)=beta*(lambda^2-1)/p,D(f)(-lambda)=p},numeric,parameters=[p],output=listprocedure);
F:=subs(init,f(x));
q:=proc(pp)
   F(parameters=[pp]);  
   F(1);
end proc;
#Testing q
q(.15);
plot(q,-.15..0.15,-1..1);
plot(q,.023..0.024,-1..1);
plot(q,.11..0.12,-1..1);
#The plots reveal approximately the values of p which make f(1) zero.

q(.0236);
q(-.0236);
F(parameters=[-.0236]);
F(1);
plot(F,-lambda..1);
F(parameters=[-.116]);
plot(F,-lambda..1);
F(1);

#An illustration

Q:=proc(pp)
   F(parameters=[pp]);  
   plot(F,-lambda..1)
end proc;
plots:-animate(Q,[p],p=-0.15..0,frames=50);

I agree that D(h) ought not return 0. Returning unevaluated would be OK.

restart;
#showstat(D);
f:=k->(y->y^k);
h:=z->f(2)(z);

# In the help page for D it is described as doing something like this:
D(h) = (q -> eval(diff(h(t),t), t=q));
#However, that is only meant to convey the idea as you see in the following (and by looking at the code for D using showstat(D)).
infolevel[D]:=5:
D(h)(w) = (q -> eval(diff(h(t),t), t=q))(w);
D: Applying  D to  h
D: Either a name or a procedure
D: Looking for D/h
D: Looking for diff/h
D: No pre-defined rules
D: Applying  D to  f
D: Either a name or a procedure
D: Looking for D/f
D: Looking for diff/f
D: No pre-defined rules
D: Applying  D to  f(2)
                            0 = 2 w
#If you redefine f in the following way (reminiscent of the way procedures returning procedures had to be handled many versions ago) then D(h) returns unevaluated
f:=proc(k) subs(_k=k,y->y^_k) end proc;

D(h)(w) = (q -> eval(diff(h(t),t), t=q))(w);
D: Applying  D to  h
D: Either a name or a procedure
D: Looking for D/h
D: Looking for diff/h
D: No pre-defined rules
D: Applying  D to  f
D: Either a name or a procedure
D: Looking for D/f
D: Looking for diff/f
D: No pre-defined rules
D: Applying  D to  subs
D: Either a name or a procedure
D: Looking for D/subs
D: Looking for diff/subs
D: No pre-defined rules
                         D(h)(w) = 2 w



You must use a list to break an existing order. Once a particular order like {a,c,b} is chosen, then the command {a,b,c}; results in {a,c,b}.
Suppose that res is the result of your fsolve command.

res:=fsolve(.....);
U:=[u1,u2,u3,u4,u5];
zip(`=`,U,eval(U,res));

# In versions later than Maple 12 you can use elementwise equality.
U=~eval(U,res);


You don't tell us what xi_M, anteil, or ableitung are, so I just made a very simplified version. With the simplified version deltaT remains a Matrix.

deltaT :=Matrix(1..19,1..6,fill=1);
i := 0;
j := 0;
for a from 1 by 1 to 19 do i:=i+1; j:=0;
for b from 1 by 1 to 6 do  j:=j+1;
deltaT(i,j):=i+j;
end do:
end do:

whattype(deltaT);

or why not just

for i from 1 to 19 do
   for j from 1 to 6 do
      deltaT(i,j):=i+j;
   end do:
end do:

whattype(deltaT);

There were typos in H, which is not a procedure/function, but just an expression containing tau.

So try this:

fQ1,fQ2,fP1,fP2 := op(subs(p,[Q1(tau),Q2(tau),P1(tau),P2(tau)])):
H := (1/8)*(fP1^2+fP2^2)/(fQ1^2+fQ2^2)+1/2*(fP1*fQ2-fQ1*fP2)+q*(fQ1^2+fQ2^2)/(1+q)-1/((1+q)*(fQ1^2+fQ2^2))-q/((1+q)*sqrt((fQ1^2-fQ2^2-1)^2+4*fQ1^2*fQ2^2))-q^2/(2*(1+q)^2):
 plot(H, tau = 0.3e-1 .. T);

First 134 135 136 137 138 139 140 Last Page 136 of 160