## 19 Badges

13 years, 105 days

## You define an ODE system...

with the independent variable eta, the dependent variables  f(eta), g(eta), phi(eta), and theta(eta), togeether with the parameters M, N, Pr, S, Sc, d, sigma, C[T], N[b], N[t], Re[r],beta^%H, and d^%H.

You seem to want plots of "Streamlines, Isotherms and Microrotations".

How are these "Streamlines, Isotherms and Microrotations". defined in terms of the variables and parameters of the ODE system listed above?

According to the graphic you supply, the "Streamlines, Isotherms and Microrotations" are parameterised by Re , Pr,  Gr and Ha. However Re, Gr, and Ha occur nowhere in the definition of the ODE system. What are these parameters?

## Since everyting...

in your worksheet appears to execute exactly as expected (see the attachecd), you are going to have to do a much better job of explaining your what your problem is

 > restart; d1 := 0.05; d2 := 0.3; AA := 0.2; BB := 0.1; PDE1 := diff(u(x, t), t) = d1*diff(u(x, t), x, x) + w(x, t)*exp(AA*u(x, t) - BB*v(x, t)); PDE2 := diff(v(x, t), t) = d2*diff(v(x, t), x, x) - w(x, t)*exp(AA*u(x, t) - BB*v(x, t)); PDE3 := 0.0001*diff(w(x, t), t) = diff(w(x, t), x) - 0.8*x + 3.3; IBC1 := u(0, t) = 1, u(1, t) = 0, u(x, 0) = piecewise(x < 0.35, -(4*x)*x + 1, 0.35 < x and x < 0.65, 1.32958 - 1.29167*x, 0.65 < x, 4*(x - 1)^2); IBC2 := v(0, t) = 0, v(1, t) = 1, v(x, 0) = piecewise(x < 0.35, (4*x)*x + 1, 0.35 < x and x < 0.65, 1.32958 - 1.29167*x, 0.65 < x, -4*(x - 1)^2); IBC3 := w(0, t) = 0.5, w(x, 0) = 1 - (0.3*x)*x; pds := pdsolve([PDE1, PDE2, PDE3], [IBC1, IBC2, IBC3], numeric, time = t, range = 0 .. 1); p1 := pds:-plot(t = 0, numpoints = 50); p2 := pds:-plot(t = 1/8, numpoints = 50, color = blue); p3 := pds:-plot(t = 1/4, numpoints = 50, color = green);
 >

Download pdeStuff.mw

## So...

something like the attached - maybe?

 > L:= [ ["O",  3.85090000,  0.45160000,  0.00120000],         ["O", -2.59990000,  1.40410000, -0.00180000],         ["N", -1.57050000, -0.71710000,  0.00010000],         ["C", -0.20660000, -0.42310000, -0.00020000],         ["C",  0.22050000,  0.90470000,  0.00040000],         ["C",  0.72980000, -1.45700000, -0.00070000],         ["C",  1.58410000,  1.19860000,  0.00020000],         ["C",  2.09330000, -1.16290000, -0.00070000],         ["C",  2.52040000,  0.16480000, -0.00030000],         ["C", -2.64850000,  0.17820000,  0.00090000],         ["C", -3.97350000, -0.54200000,  0.00100000],         ["H", -0.44360000,  1.75770000,  0.00120000],         ["H",  0.41130000, -2.49630000, -0.00100000],         ["H", -1.80100000, -1.70860000,  0.00010000],         ["H",  1.90530000,  2.23700000,  0.00090000],         ["H",  2.81800000, -1.97260000, -0.00080000],         ["H", -4.06550000, -1.14630000, -0.90580000],         ["H", -4.79040000,  0.18440000,  0.02880000],         ["H", -4.04450000, -1.18860000,  0.88020000],         ["H",  3.96500000,  1.41760000,  0.00170000]       ]:   a1:={seq( `if`( evalb(L[j,1]="H"), j, NULL), j=1..numelems(L))};   S:={ {1,9}, {1,20},  {2,10},   {3,4},  {3,10},   {3,14}, {4,5},        {4,6}, {5, 7},  {5,12},   {6,8},  {6,13},   {7,9} ,{7,15},        {8,9}, {8,16}, {10,11}, {11,17}, {11,18}, {11,19}      }:   a2:={seq(`if`(`intersect`(j, a1)={}, j, NULL), j in S)};
 (1)
 >

Download manip.mw

## Change...

the drawGraph command to

`DrawGraph(G, style=spring, showweights=true);`

This works for me, although you may have to start playing with font sizes in order to make things legible

My reading of the relevant section of the help page reproduced (below)

showweights=truefalse
The display of edge weights can be forced or suppressed with the option showweights=true and showweights=false, respectively. By default, the edge weights are displayed when G is a weighted graph with fewer than 46 edges.

is that if the graph has more than 46 edges the display of weights will be suppressed, although interestingly your graph *only* has 39 edges

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

Download FT.mw

## 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]`

## What advantages...

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

Download ODEnotPDE3.mw

## Read my previous responses...

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.

Please remeber that your original code was so full of logical and syntax errors that I had to guess your intention. It is quite likely that I made a bad guess :-(

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

Download ODEnotPDE2.mw

## 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

 2 3 4 5 6 7 8 Last Page 4 of 206
﻿