Preben Alsholm

13728 Reputation

22 Badges

20 years, 250 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Notice that by the assignments to x[i+1/2] and x[i+1] you are (implicitly) creating a table with name x. Try
eval(x);

Notice also that those assignments don't help any if  i  is concrete: Try
x[1+1/2];
i:=7;
x[i+1]-x[i];

There is the same "problem" with phi. Not a real problem, but you mentioned aesthetics.
The code does become a problem if i is in fact assigned a value (as I did above in i:=7;) say just before the assignments to the K's.

Notice that by the assignments to x[i+1/2] and x[i+1] you are (implicitly) creating a table with name x. Try
eval(x);

Notice also that those assignments don't help any if  i  is concrete: Try
x[1+1/2];
i:=7;
x[i+1]-x[i];

There is the same "problem" with phi. Not a real problem, but you mentioned aesthetics.
The code does become a problem if i is in fact assigned a value (as I did above in i:=7;) say just before the assignments to the K's.

restart;
interface(rtablesize=infinity):
N:=10:
K:=Matrix(N-1):
for i from 1 to N do x[i]:=x[i-1]+h end do:
for i from 1 to N do
 phi[i]:=piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]end do:
for i from 2 to N-2 do
    K[i,i]:=int(diff(phi[i], t)*diff(phi[i], t), t=x[i]..x[i+1]) assuming h>0;
    K[i,i+1]:=int(diff(phi[i], t)*diff(phi[i+1], t), t=x[0]..x[N]) assuming h>0;
    K[i-1,i]:=K[i,i+1];
 end do:
K;

Since I don't know what your intention is I cannot say if the statement K[i-1,i]:=K[i,i+1]; is correct or not.

There is an obvious pattern, so maybe you need not do any other calculations than

f:=piecewise(t<=a-h,0,t<a, (t-a)/h, t<=a+h, (a+h-t)/h, 0);
int(diff(f, t)*diff(f, t), t=a..a+h) assuming h>0;
int(diff(f, t)*diff(eval(f,a=a+h), t), t=a..a+h) assuming h>0;
#The same as:
int(diff(f, t)*diff(eval(f,a=a+h), t), t=a-h..a+2*h) assuming h>0;

restart;
interface(rtablesize=infinity):
N:=10:
K:=Matrix(N-1):
for i from 1 to N do x[i]:=x[i-1]+h end do:
for i from 1 to N do
 phi[i]:=piecewise(t>=x[i-1] and x[i]>t, (t-x[i])/(h), x[i]end do:
for i from 2 to N-2 do
    K[i,i]:=int(diff(phi[i], t)*diff(phi[i], t), t=x[i]..x[i+1]) assuming h>0;
    K[i,i+1]:=int(diff(phi[i], t)*diff(phi[i+1], t), t=x[0]..x[N]) assuming h>0;
    K[i-1,i]:=K[i,i+1];
 end do:
K;

Since I don't know what your intention is I cannot say if the statement K[i-1,i]:=K[i,i+1]; is correct or not.

There is an obvious pattern, so maybe you need not do any other calculations than

f:=piecewise(t<=a-h,0,t<a, (t-a)/h, t<=a+h, (a+h-t)/h, 0);
int(diff(f, t)*diff(f, t), t=a..a+h) assuming h>0;
int(diff(f, t)*diff(eval(f,a=a+h), t), t=a..a+h) assuming h>0;
#The same as:
int(diff(f, t)*diff(eval(f,a=a+h), t), t=a-h..a+2*h) assuming h>0;

You have the very same issue now on the x-axis.
Try

loglogplot(j(x), x = 1/1000 .. 7*10^7, thickness = 1);

You have the very same issue now on the x-axis.
Try

loglogplot(j(x), x = 1/1000 .. 7*10^7, thickness = 1);

As always you are so much more likely to get a useful response if you provide details, in this case the equations as text or as an uploaded worksheet.

Obviously you didn't get plotting Sin[x*Pi]*Exp[A*x] to work for any value of A if you used that syntax.
So which command did you actually ask Maple to execute?

@arunousephthoma My concern was whether you have the correct branches of the different square roots. So the basic question is: How did you end up with square roots? The original problem most likely didn't contain any?

And by the way, you may try

plots:-complexplot3d(m*10^(-24),x =1*10^5+1*10^4*I..1*10^7+1*10^7*I);

What is the original problem? Obviously this didn't come out of the blue.
Which symbolic transformations took place before you ended up with this?
The appearance of square roots of complex numbers suggests the question: Are you sure to have the right branch of the square root?
The square root in Maple is the principal one.

@navid Are there any particular values for the parameters for which this happens? Could you give an example.

Like e.g.

param:={n=50,lambda=.1,p=1/3};

for which I got 50 roots.

@trottaleonardo OK, so you tried a different approach. But your code makes a 1 x nops(Z) matrix.
Revise it to (assuming w is the same as your Z):

convert(w,list);
`[]`~(%); #To turn the list into a listlist
convert(%,Matrix);

In the example you gave w was a sum of terms (type `+`). I'm assuming this is also the case with Z. If nops(Z) = 2 and Z is of type `+` then there are only 2 terms. So you get only two lines.

By the way, did you actually try my first code on your own example?

@trottaleonardo OK, so you tried a different approach. But your code makes a 1 x nops(Z) matrix.
Revise it to (assuming w is the same as your Z):

convert(w,list);
`[]`~(%); #To turn the list into a listlist
convert(%,Matrix);

In the example you gave w was a sum of terms (type `+`). I'm assuming this is also the case with Z. If nops(Z) = 2 and Z is of type `+` then there are only 2 terms. So you get only two lines.

By the way, did you actually try my first code on your own example?

Your code exactly as written (except for supplying a missing semicolon at the end) worked for me.  I'm using ghostscript and ghostview for viewing. 

Try

Digits:=25;

and a much smaller range (0..0.0001).

Added later:

In order to see that X is not constant (although it appears to be at first sight) you can plot like this:

rg:=0..0.0001;
X0:=10^20:
init:=R(0)=0,X(0)=X0:
res:=dsolve({sys,init},[X(x),R(x)],numeric,range=rg);
plots:-odeplot(res,[x,R(x)],rg);
plots:-odeplot(res,[x,ln(X(x)/X0)],rg,thickness=4,view=-1.5e-17..0);

#Alternatively, make a simple change of variableX(x) = X0 + xi(x):

X0:='X0':
sys2:=PDEtools:-dchange({X(x)=X0+xi(x)},{sys},[xi]);
Digits:=25:
X0:=10^20:
init2:=xi(0)=0,R(0)=0:
rg:=0..0.0001;
res:=dsolve(sys2 union {init2},[xi(x),R(x)],numeric,range=rg);
plots:-odeplot(res,[x,R(x)],rg);
plots:-odeplot(res,[x,xi(x)],rg,thickness=4);



First 181 182 183 184 185 186 187 Last Page 183 of 230