Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

My preference would be

evalindets(sol,specfunc(anything,{sin,cos}),s->op(0,s)(omega*t));

But your way could be made similarly:

evalindets(sol,specfunc(anything,{sin,cos}),x->applyrule(eq,x));

@mapleq2013 x(t)-1 is a rootfinding trigger, meaning that when during integration x(t) - 1 hits zero the action (in this case 'halt') is activated.
So in the example there are two events both of the type simple rootfinding and both having action 'halt'.

Here is a simple example of a conditional trigger:

sys2:=diff(x(t),t)= y(t) , diff(y(t),t)= -x(t);  
res2:=dsolve({sys2,x(0)=1/2,y(0)=0},numeric,events=[[[x(t),y(t)>0],halt]]);
plots:-odeplot(res2,[t,x(t)],0..5);

In this one integration stops when x(t) = 0 provided that also y(t) > 0.

@mapleq2013 x(t)-1 is a rootfinding trigger, meaning that when during integration x(t) - 1 hits zero the action (in this case 'halt') is activated.
So in the example there are two events both of the type simple rootfinding and both having action 'halt'.

Here is a simple example of a conditional trigger:

sys2:=diff(x(t),t)= y(t) , diff(y(t),t)= -x(t);  
res2:=dsolve({sys2,x(0)=1/2,y(0)=0},numeric,events=[[[x(t),y(t)>0],halt]]);
plots:-odeplot(res2,[t,x(t)],0..5);

In this one integration stops when x(t) = 0 provided that also y(t) > 0.

@mtj4u The result from the integration over the finite interval r..infinity is the result of the integral from 0 to infinity of f with an error less than the value of the sum of epsilon and the integral of g from r to infinity. In the code epsilon = 1e-10 and really means relative error, but since the result is close to 1, there is no difference between relative and absolute error. 

Since I observed some error messages in your worksheet I upload my version with r = 10^6.

MaplePrimesIntegral1.mw

Edited: 'sum of' has replaced 'maximum of'. The text in the worksheet has been adjusted accordingly.

NB. You are using Maple 14 it appears. I'm using Maple 16. The maxintervals option for the methods _d01ajc and _d01akc is not available in Maple 12 I just found out. I have no access to Maple 14, but that option is available in Maple 15.

@mtj4u The result from the integration over the finite interval r..infinity is the result of the integral from 0 to infinity of f with an error less than the value of the sum of epsilon and the integral of g from r to infinity. In the code epsilon = 1e-10 and really means relative error, but since the result is close to 1, there is no difference between relative and absolute error. 

Since I observed some error messages in your worksheet I upload my version with r = 10^6.

MaplePrimesIntegral1.mw

Edited: 'sum of' has replaced 'maximum of'. The text in the worksheet has been adjusted accordingly.

NB. You are using Maple 14 it appears. I'm using Maple 16. The maxintervals option for the methods _d01ajc and _d01akc is not available in Maple 12 I just found out. I have no access to Maple 14, but that option is available in Maple 15.

@PlpPlp I changed things a little to reflect my understanding of the problem. I use worksheet interface and 1D input (so much easier for me).

DeltaF := K[1]*m[1]^(2)*m[2]^(2)+(K[1]+(B[1]^(2))/(2* c[11])-(B[2]^(2))/(2 *c[44]))*(m[1]^(2)+m[2]^(2))*m[3]^(2)+K[2]*m[1]^(2)*m[2]^(2)*m[3]^(2)+B[1]*(u[m1]*m[1]^(2)+u[m2]*m[2]^(2))+B[2]*u[m6]*m[1]^(2)*m[2]^(2)-B[1]*((B[1])/(6 *c[11])+(c[12])/(c[11])*(u[m1]+u[m2]))*m[3]^(2)+1/(2)*mu[0]*M[s]^(2)*(N[1]*m[1]^(2)+N[2]*m[2]^(2)+N[3]*m[3]^(2))-mu[0]*M[s]*(H[1]*m[1]+H[2]*m[2]+H[3]*m[3]);

H[eff1] := -diff(DeltaF, m[1])/mu[0];
H[eff2] := -diff(DeltaF, m[2])/mu[0];
H[eff3] := -diff(DeltaF, m[3])/mu[0];

expr := collect(mu[0]*(m[2]*H[eff3]-m[3]*H[eff2]),[m[1],m[2],m[3]],distributed,factor);

#Now linearising about m = m0. The number 2 at the end means that the orders neglected are 2 and above.

mtaylor(expr, [m[1]=m[1,0],m[2]=m[2,0],m[3]=m[3,0] ], 2);

@PlpPlp I changed things a little to reflect my understanding of the problem. I use worksheet interface and 1D input (so much easier for me).

DeltaF := K[1]*m[1]^(2)*m[2]^(2)+(K[1]+(B[1]^(2))/(2* c[11])-(B[2]^(2))/(2 *c[44]))*(m[1]^(2)+m[2]^(2))*m[3]^(2)+K[2]*m[1]^(2)*m[2]^(2)*m[3]^(2)+B[1]*(u[m1]*m[1]^(2)+u[m2]*m[2]^(2))+B[2]*u[m6]*m[1]^(2)*m[2]^(2)-B[1]*((B[1])/(6 *c[11])+(c[12])/(c[11])*(u[m1]+u[m2]))*m[3]^(2)+1/(2)*mu[0]*M[s]^(2)*(N[1]*m[1]^(2)+N[2]*m[2]^(2)+N[3]*m[3]^(2))-mu[0]*M[s]*(H[1]*m[1]+H[2]*m[2]+H[3]*m[3]);

H[eff1] := -diff(DeltaF, m[1])/mu[0];
H[eff2] := -diff(DeltaF, m[2])/mu[0];
H[eff3] := -diff(DeltaF, m[3])/mu[0];

expr := collect(mu[0]*(m[2]*H[eff3]-m[3]*H[eff2]),[m[1],m[2],m[3]],distributed,factor);

#Now linearising about m = m0. The number 2 at the end means that the orders neglected are 2 and above.

mtaylor(expr, [m[1]=m[1,0],m[2]=m[2,0],m[3]=m[3,0] ], 2);

You have posted several of those by now. However, I fail to see the connection to mathematical software like Maple.

Maybe uploading the worksheet with your example could help make the problem clear.

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);

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