## 13776 Reputation

14 years, 62 days

## Maybe...

but the following paragraph (emphasis added) from https://en.wikipedia.org/wiki/Autocorrelation appeared significant.

Different fields of study define autocorrelation differently, and not all of these definitions are equivalent. In some fields, the term is used interchangeably with autocovariance.

## Hmmmmm...

@Carl Love
Some excerpts from Maple help (emphasis added occasionally)

A set is an unordered sequence of distinct expressions enclosed in braces {}, representing a set in the mathematical sense.

Note:  This order could be changed in future versions of Maple, and/or objects could be represented in a different way in future versions. You should not assume that sets will be maintained in any particular order.

Please  feel free to use "order"  with sets!

## When using the geometry() package...

There is a bit of a trade-off to be considered. If everyy input is expressed in exact numbers ( eg rationals or integers) then all calculations will be performed using exact numbers. Depending on the complexity of the calculations you are performing, this can lead to very complicated expressions, which *may* skiw down processing. One can circumvent this problem by supplying iall numerica inputs as floats (see the attached, which will execute "faster" than my original response).

There are however drawbacks to using floats. Certain commands within the geometry package require computing exact equality between two different expressions (consider for example, IsOnCircle() or IsOnLine() ). As a general programming rule, testing for equality between two floating-point numbers, computed in two different calculations is a big no-no. Being off by 1 in the 10th decimal place will kill you.

You have a choice between time and code robustness - pick one because you can't have both!!

## Things I don't understand...

@arashghgood

1. Your ode system can be solved analytically (see the attached), so why are you doing it numerically???
2. Your worksheet states "u_hat[i]s are in Fourier space. I have to apply inverse Fourier transform and plot the u(t)". But your worksheet has the independent variable of all u_hat[i] as 't'. For clarity, should this be u_hat[i](omega) ? I know this may just be a notation issue, but you should make it really clear whether the uhat[i] are frequency- or time-domain functions. That way it is clear whether one needs to use a Fourier transfrom or an inverse Fourier transform to generate the correspondin time- or frequency-domain functions.

 > restart;
 > with(DEtools):
 > # define the complex numbers list1 and list2 list1 := [ [0., -5.496799068*10^(-15)-0.*I], [.1, 5.201897725*10^(-16)-1.188994754*I], [.2, 6.924043163*10^(-17)-4.747763855*I], [.3, 2.297497722*10^(-17)-10.66272177*I], [.4, 1.159126178*10^(-17)-18.96299588*I] ] ; list2 :=[ [0., -8.634351786*10^(-7)-67.81404036*I], [.1, -0.7387644021e-5-67.76491234*I], [.2, -0.1433025271e-4-67.59922295*I], [.3, -0.2231598645e-4-67.25152449*I], [.4, -0.3280855430e-4-66.56357035*I] ];
 (1)
 > # define the differential equation system ode_system := []; ic_system :=[]; for i from 1 to 5 do   ode_system := [op(ode_system), diff(u_hat[i](t),t,t) + I*(list1[i][2]+list2[i][2])*diff(u_hat[i](t),t) - list1[i][2]*list2[i][2]*u_hat[i](t) = 0];   ic_system := [op(ic_system), u_hat[i](0) = 1, D(u_hat[i])(0) = 0] end do:
 (2)
 > #sol := dsolve( map(op, [ode_system, ic_system]), numeric, output = listprocedure );
 (3)
 > #If it is correct, u_hat[i]s are in Fourier space. I have to apply inverse Fourier transform and plot the u(t). #Is the dsolve above correct? #How do next steps?
 > # # Why is OP doing this problem numerically? #   sol:=evalf~(dsolve( map(op, [ode_system, ic_system]))):   eval( u_hat[1](t), sol);   eval( u_hat[2](t), sol);   eval( u_hat[3](t), sol);   eval( u_hat[4](t), sol);   eval( u_hat[5](t), sol);
 (4)
 >
 >

## Post...

the complete worksheet using the big green up-arrow in the Mapleprimes toolbar

## numeric...

`Error, (in dsolve) too many arguments; some or all of the following are wrong: [{u[1](t), u[2](t), u[3](t), u[4](t), u[5](t)}, numerical, output = listprocedure]`

should be numeric, as in

` [{u[1](t), u[2](t), u[3](t), u[4](t), u[5](t)}, numeric, output = listprocedure]`

does this code offer compared with the NyQuistPlot() command in Maple's DynamicSystems() package?

See the help at ?NyQuistPlot

## The latest code you present...

does not execute in Maple 2022. Too many errors for me to debug

n my Maple 2018, some of the execution groups do not produce the output which you display. I have no idea how you achieved this - and I am no longer prepared to guess.

The only thing I can suggest is that you stop writing code, and instead analyse the attached very carefully. When you understand it, you *may* be able to point out whether or not I have made a logical error.

It might also be useful to supply the original ODEs so that I  could run Maple's numerical ODE solver, to compare results with this implementation of the "hpm method"

 > restart:   N:= 3: M:= 0.1: Kp:= 0.1: Gr:= 0.01: Gc:= 0.01:   Pr:= 1: S:= 0.01: Sc:= 0.78: Kc:= 0.01: La:= 1:   f__N:=  add((p^(i))*f[i](x), i = 0 .. N):   Theta__N:=  add((p^(i))*Theta[i](x), i = 0 .. N):   Phi__N:= add((p^(i))*Phi [i] (x), i = 0 .. N):   odeSys:= { (1-p)*(diff(f__N, x, x, x))+p*(diff(f__N, x, x, x)+(1/2)*(diff(f__N, x, x))*f__N-(M^2+Kp)*(diff(f__N, x)-La)+Gr*Theta__N+Gc*Phi__N),              (1-p)*(diff(Theta__N, x, x))/Pr+p*((diff(Theta__N, x, x))/Pr+(1/2)*(diff(Theta__N, x))*f__N+S*Theta__N),              (1-p)*(diff(Phi__N, x, x))/Sc+p*((diff(Phi__N, x, x))/Sc+(1/2)*(diff(Phi__N, x))*f__N+Kc*Phi__N)            }:   cond:= f[0](0) = 0,     D(f[0])(0) = 0,  D(f[0])(5) = 1,          Theta[0](0) = 1, Theta[0](5) = 0,          Phi[0](0) = 1,   Phi[0](5) = 0:   ans:={}:   for k from 0 by 1 to N do       ans:= `union`             ( ans,               dsolve               ( { eval                   ( coeff~(odeSys, p, k),                     ans                   )[],                   cond                 }               )            ):       cond:= f[k+1](0) = 0, D(f[k+1])(0) = 0, (D(f[k+1]))(5) = 0,              Theta[k+1](0) = 0, Theta[k+1](5) = 0,              Phi[k+1](0) = 0, Phi[k+1](5) = 0:   od:   getfunc:= z-> z=eval( z, evalf~(ans)):
 > for k from 0 to N do       getfunc( f[k](x));   od;   for k from 0 to N do       getfunc( Theta[k](x));   od;   for k from 0 to N do       getfunc( Phi[k](x));   od;
 (1)
 > interface(rtablesize=25):   Matrix( [ [ eta, f(eta), Theta(eta), Phi(eta) ],               seq               ( [ j,                   eval(add( rhs(getfunc(f[k](x))), k=0..N), x=j),                   eval(add( rhs(getfunc(Theta[k](x))), k=0..N), x=j),                   eval(add( rhs(getfunc(Phi[k](x))), k=0..N), x=j)                 ],                 j=0..5, 0.25               )           ]       );   interface(rtablesize=25):
 (2)
 >

for example

If you want the actual data at a finer resolution, then just change the step size(s) in the final execution group.

and

just change the endpoints of the seq commands for example
Matrix([seq([seq(interfunc(i,j),i=0..1,0.001)],j=1..0, -0.001)]);

## As close as I can get...

to the format you require (which is barely legible!). A quich check suggest that successive the equations in the attached  for 'Theta' and 'Phi' agree with you 'picture' , but those in 'f' do not. You should check this carefully.

 > restart:   N := 4: M := .1: Kp := .1: Gr := 0.1e-1: Gc := 0.1e-1: Pr := 1: S := 0.1e-1: Sc := .78: Kc := 0.1e-1: La := 1:   f__N:=  sum((p^(i))*f[i](x), i = 0 .. N) :   Theta__N:=  sum((p^(i))*Theta[i](x), i = 0 .. N) :   Phi__N:= sum((p^(i))*Phi [i] (x), i = 0 .. N):   HPMEq1 := (1-p)*(diff(f__N, x, x, x))+p*(diff(f__N, x, x, x)+(1/2)*(diff(f__N, x, x))*f__N-(M^2+Kp)*(diff(f__N, x)-La)+Gr*Theta__N+Gc*Phi__N):   HPMEq2 := (1-p)*(diff(Theta__N, x, x))/Pr+p*((diff(Theta__N, x, x))/Pr+(1/2)*(diff(Theta__N, x))*f__N+S*Theta__N):   HPMEq3 := (1-p)*(diff(Phi__N, x, x))/Sc+p*((diff(Phi__N, x, x))/Sc+(1/2)*(diff(Phi__N, x))*f__N+Kc*Phi__N):   for i from 0 to N do       equ[1][i] := coeff(HPMEq1, p, i) = 0:   end do:   for i from 0 to N do       equ[2][i] := coeff(HPMEq2, p, i) = 0:   end do:   for i from 0 to N do       equ[3][i] := coeff(HPMEq3, p, i) = 0:   end do:   cond[1][0] := f[0](0) = 0, (D(f[0]))(0) = 0, Theta[0](0) = 1, Phi[0](0) = 1, Theta[0](5) = 0, Phi[0](5) = 0, (D(f[0]))(5) = 1:   for j to N do        cond[1][j] := f[j](0) = 0, (D(f[j]))(0) = 0, Theta[j](0) = 0, Phi[j](0) = 0, Theta[j](5) = 0, Phi[j](5) = 0, (D(f[j]))(5) = 0:   end do:   ans:=dsolve({equ[1][0], equ[2][0], equ[3][0], cond[1][0]}):   for k from 1 by 1 to 4 do       ans:=`union`( ans, dsolve( { eval({equ[1][k], equ[2][k], equ[3][k]}, ans)[], cond[1][k]})):   od:   for j from 0 by 1 to 4 do       f[j](x)=eval( f[j](x), evalf~(ans));   od;   for j from 0 by 1 to 4 do       Theta[j](x)=eval( Theta[j](x), evalf~(ans));   od;   for j from 0 by 1 to 4 do       Phi[j](x)=eval( Phi[j](x), evalf~(ans));   od;
 (1)
 >

## To change the range...

just change the endpoints of the seq commands for example

Matrix([seq([seq(interfunc(i,j),i=0..1,0.001)],j=1..0, -0.001)]);

will generate a 1000x1000 matrix of values from 0..1 in both directions, in steps of 0.001

## There is no such thing...

as a "dat file"

If you want the actual data at a finer resolution, then just change the step size(s) in the final execution group. You can then transfer the reults matrix anywhere you want using ExportMatrix()

## So...

Post the worksheet which produces the error (Use big green uparrow in the Mapleprimes toolbar)

## Cannot reproduce...

The Matlab code

```>> M=rand(1800,1);
>> save('C:\Users\TomLeslie\Desktop\test.mat', 'M')```

creates an 1800x1 matrix on my desktop.

Both of the maple commands shown below read it

```ImportMatrix("C:/Users/TomLeslie/Desktop/test.mat", source=Matlab, output=matrices);
ImportMatrix("C:/Users/TomLeslie/Desktop/test.mat", source=MATLAB, output=matrices);
```

I suppose it *might* be a version issue.  Which version Matlab and Maple are you running?

(I'm using Maple 2022.2 and Matlab R2022b)

## Hmmm...

fdiff() seems to be doing something pretty crazy - and I don't understand why, so, in the attached, I have changed the method of computing the numerical derivative - the answer looks a lot more plausible

At the moment I can't upload a worksheet here, because someone has broken the big green up-arrow! However you need to replace the code

```Tv := eval(T(y, t), (pds[i]):-value(t = 1, output = listprocedure));
f := yv -> (1 + Nr*(Tv(yv) + 1))*fdiff(Tv(y), y = yv, workprec = 1.5);
plt1[i] := plot('f(y)', y = 0 .. 1, color = AA[i], size = [1000, 600], axes = boxed, axesfont = ["Arial", 14, Bold], thickness = 3);```

with

```h := 0.001;
Tv := eval(T(y, t), (pds[i]):-value(t = 1, output = listprocedure));
f := yv -> (1 + Nr*(Tv(yv) + 1))*(Tv(yv + h) - Tv(yv))/h;
plt1[i] := plot('f(y)', y = 0 .. 1, color = AA[i], size = [1000, 600], axes = boxed, axesfont = ["Arial", 14, Bold], thickness = 3);```
 4 5 6 7 8 9 10 Last Page 6 of 207
﻿