Preben Alsholm

13728 Reputation

22 Badges

20 years, 250 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Here is a version which I believe is sound. It adds and subtracts series expansions of sin and cos to make the individual integrals convergent. It is essential that the sum of the terms in L1 be 0. That follows from knowing that the original integrand has no singularity at zero.

The code should work for even as well as odd values of N.
It is still slow for all but one digit values of N.

restart;
with(IntegrationTools):
A:=n->Int(mul(sin(x/2^k)*2^k/x, k = 0 .. n), x = 0 .. infinity);
#N:=9:
N:=4:
Change(A(N),x=2^N*t):
B:=combine(%):
f1:=GetIntegrand(B):
series(numer(f1),t=0,N+1); #terms of order up to and including N are zero
L:=convert(expand(f1,sin,cos),list):
sn:=unapply(convert(series(sin(x),x=0,N+1),polynom),x):
cs:=unapply(convert(series(cos(x),x=0,N), polynom),x):
L1:=eval(L,{sin = sn,cos=cs}):
simplify(convert(L1,`+`)); #Will return 0.
map(int,L-L1,t=0..infinity):
convert(%,`+`);

Here is a version which I believe is sound. It adds and subtracts series expansions of sin and cos to make the individual integrals convergent. It is essential that the sum of the terms in L1 be 0. That follows from knowing that the original integrand has no singularity at zero.

The code should work for even as well as odd values of N.
It is still slow for all but one digit values of N.

restart;
with(IntegrationTools):
A:=n->Int(mul(sin(x/2^k)*2^k/x, k = 0 .. n), x = 0 .. infinity);
#N:=9:
N:=4:
Change(A(N),x=2^N*t):
B:=combine(%):
f1:=GetIntegrand(B):
series(numer(f1),t=0,N+1); #terms of order up to and including N are zero
L:=convert(expand(f1,sin,cos),list):
sn:=unapply(convert(series(sin(x),x=0,N+1),polynom),x):
cs:=unapply(convert(series(cos(x),x=0,N), polynom),x):
L1:=eval(L,{sin = sn,cos=cs}):
simplify(convert(L1,`+`)); #Will return 0.
map(int,L-L1,t=0..infinity):
convert(%,`+`);

@Markiyan Hirnyk I overlooked the fact that the sine integrals in the first part of my answer are also divergent.
I have inserted notes in my answer.
I prefer correct reasoning with a wrong result over wrong reasoning with a correct result.

@Markiyan Hirnyk I overlooked the fact that the sine integrals in the first part of my answer are also divergent.
I have inserted notes in my answer.
I prefer correct reasoning with a wrong result over wrong reasoning with a correct result.

Have you managed to make this work in a simpler case?
For instance with r0 =1  and rf = 10 and number of partitions = 2?

@reemeaaaah In order to have total control over the function calls I thought it safer to use a bare bone version of Euler's method.
I made two versions. Euler is the usual and requires input of the function f in y' = f(t,y).
Euler2 requires input of a function f of 3 variables y' = f(t,y,r), where the third argument is meant for the random input.
Your particular normal distribution is hardcoded into Euler2.
For what it is worth, here it is:

StochasticEuler.mw

@reemeaaaah In order to have total control over the function calls I thought it safer to use a bare bone version of Euler's method.
I made two versions. Euler is the usual and requires input of the function f in y' = f(t,y).
Euler2 requires input of a function f of 3 variables y' = f(t,y,r), where the third argument is meant for the random input.
Your particular normal distribution is hardcoded into Euler2.
For what it is worth, here it is:

StochasticEuler.mw

@reemeaaaah I have a habit of editing my answers a little too often. I just added some lines.
My knowledge of statistics is basically confined to the ability of spelling the word, in two languages, though.

@reemeaaaah I have a habit of editing my answers a little too often. I just added some lines.
My knowledge of statistics is basically confined to the ability of spelling the word, in two languages, though.

What effect do the extra arguments 100 and confidence=0.95 in the assignment
R := RandomTools:-Generate(
      distribution(Normal(-0.2e-1, 0.4e-1)),100,confidence=0.95, makeproc = true);
have?
It seems to me that they are just ignored. Try e.g. misspelling them or changing them in another way.
Before and after try
R();

The following syntax generates a list of 10 numbers
R := RandomTools:-Generate(list(distribution(Normal(-0.2e-1, 0.4e-1)),10));

and this
R := RandomTools:-Generate(list(distribution(Normal(-0.2e-1, 0.4e-1)),10),makeproc=true);
makes a procedure which at each invocation does the same thing:
R();

I was writing in frustration after having answered and afterwards edited the answer to:

http://www.mapleprimes.com/questions/143595-E-With-One-Random-Variable

At least my frustration got vented.

I was writing in frustration after having answered and afterwards edited the answer to:

http://www.mapleprimes.com/questions/143595-E-With-One-Random-Variable

At least my frustration got vented.

@reemeaaaah The MaplePrimes editor very often leaves out everything on a line after <.
That is quite a nuissance.
I corrected it. I hope it is OK now.

@reemeaaaah The MaplePrimes editor very often leaves out everything on a line after <.
That is quite a nuissance.
I corrected it. I hope it is OK now.

Well, I did look at it.

Did you try

evalf[20](eval(yg,p=0.3456));

evalf[20](eval(wronski,p=0.3456));

There is a huge discrepancy. My guess is that your wronski is wrong. I didn't try to follow the argument.

First 179 180 181 182 183 184 185 Last Page 181 of 230