1 years, 18 days

You are right...

You are right... The theta(z,0) = thetab*sin(pi/2*z). It was a typo. The numerical analysis would serve as a background to solve the full system. At moment we are dealing with a linear theory and we are about to consider the nonlinear analysis, which is why I need the numerical scheme for this case first before thinking of the full system.

Two ways to do this...

We use theta(z,0) = pi/2 or theta(z,0) = thetab*sin(pi/2)*z if you like. The analytical analysis supposed the I.C  theta(z,0) = thetab*sin(pi/2)*z for a very small value of thetab. Intuitively, theta(z,0) = thetab*sin(pi/2)*z IC is not compatible with the theta(d,t)=0 BC, but for a very small value of thetab, especially when thetab is almost zero, it does satisfy the BC theta(d,t)=0 BC.

Analytical Solution Known...

@dharr

I have solved it analytically. But we wanna solve the same problem numerically and compare the analytical and numerical solutions.

theta and v are functions of z and t...

@nm Sorry, theta and v are functions of z and t

Mistake corrected...

@nm Hi,

Sorry my for the mistakes. I have corrected them.  Thanks

@tomleslie

Sorry, I send the wrong file to you. The f1..f8 are supposed to be dS/dt, dE/dt, dI/dt, dH/dt, dR/dt, dD/dt, dB/dt and dP/dt, where S, E, I, H, R, D, B and P are functions of t. The bifurcation parameter should be d_1. If I'm done with the analysis, I will get back to you for the bifurcation diagram.

Thank you.

Okay, thank you

Too many recursion and Float(undefined)...

I tried to plot 101 x 101 grids but received an error: too many recursion.  I then tried it on 51 x 101 grids, and it works but the I received undefine values for both the maximum and minimum. See below

> ## Assign the min and max value of pmin to minv and maxv

#(minv,maxv):=[min,max](data[gridji,4])[];

conts := [Float(undefined), Float(undefined), Float(undefined), Float(undefined), Float(undefined), Float(undefined), Float(undefined), Float(undefined), Float(undefined), Float(undefined)]

 > restart:
 > doCalc:= proc(xi, u)  #u is the \bar(H): normalize magnetic field magnitude,                           # where H = bar(H)*H__a                  option remember;                  # Import Packages                  uses ArrayTools, Student:-Calculus1, LinearAlgebra,                       ListTools, RootFinding, plots, ListTools:                  local gamma__1:= .1093,                        alpha__3:= -0.1104e-2,                        k__1:= 6*10^(-12),                        d:= 0.2e-3,                        theta0:= 0.1e-3,                        eta__1:= 0.240e-1, chi:= 1.219*10^(-6),                        alpha:= 1-alpha__3^2/(gamma__1*eta__1),                        theta_init:= theta0*sin(Pi*z/d),                        H__a:= Pi*sqrt(k__1/chi)/d,                        H := u*H__a,                             c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2)),                        w := chi*H^2*eta__1*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),                        n:= 20,                        g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,                        qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,                        soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,                        AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar;                  if not [xi,u]::[numeric,numeric] then                     return 'procname'(args);                  end if; # Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes                  g:= q-(1-alpha)*tan(q)- (w*q + c)*tan(q):                  f:= subs(q = x+I*y, g):                  b1:= evalc(Re(f)) = 0:                  b2:= evalc(Im(f)) = 0:                  qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):                  OddAsymptotes:= Vector[row]([seq(evalhf(1/2*(2*j + 1)*Pi), j = 0 .. n)]); # Compute Odd asymptote                  ModifiedOddAsym:= abs(-~(OddAsymptotes, qstar));                  qstarTemporary:= min(ModifiedOddAsym);                  indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);                  qstar2:= OddAsymptotes(indexOfqstar2); # Compute Odd asymptote                  AreThereComplexRoots:= type(true, 'truefalse');                  try                       soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});                       soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});                       qcomplex1:= subs(soln1, x+I*y);                       qcomplex2:= subs(soln2, x+I*y);                  catch:                        AreThereComplexRoots:= type(FAIL, 'truefalse');                  end try; # Compute the rest of the Roots                  #OddAsymptotes:= Vector[row]([seq(evalhf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);                  AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));                  if   AreThereComplexRoots                  then                       gg := [ qcomplex1,                               qcomplex2,                               op(Roots(g, q = 0.1e-3 .. AllAsymptotes[n], numeric))                             ];                  elif not AreThereComplexRoots                  then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[n], numeric))                            ];                  end if: # Remove the repeated roots if any & Redefine n                  qq:= MakeUnique(gg):                  m:= numelems(qq): ## Compute the time constants                  for i to m do                  p[i] := evalf(gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2)): if not p[i]::complex(numeric) then print(p[i], [xi,u], qq[i]); end if;                  end do:                 Digits := 15; ## Compute θ_n from initial conditions                 for i to m do                 Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):                 end do: ## Compute integral coefficients for off-diagonal elements θ_n matrix                 printlevel := 2:                 for i to m do                     for j from i+1 to m do                         b[i, j] := evalf(Int(Efun[i]*Efun[j], z = 0 .. d)):                         if not b[i, j]::complex(numeric) then                             b[i, j] := evalf[15](Int(Efun[i]*Efun[j], z = 0 .. d,                                                  digits=15, method=_Dexp, epsilon=1e-12)):                             if not b[i, j]::complex(numeric) then                                print("A",[Efun[i]*Efun[j], z = 0 .. d]);                             end if;                         end if;                         b[j, i] := b[i, j]:                         aa[i, j] := b[i, j]:                         aa[j, i] := b[j, i]:                     end do :                 end do: ## Calculate integral coefficients for diagonal elements in theta_n matrix                 for i to m do                    aa[i, i] := evalf(Int(Efun[i]^2, z = 0 .. d)):                    if not aa[i, i]::complex(numeric) then                       aa[i, i] := evalf[15](Int(Efun[i]^2, z = 0 .. d,                                                 digits=15, epsilon=1e-12));                       if not aa[i, i]::complex(numeric) then                          print("B",[Efun[i]^2, z = 0 .. d]);                       end if;                    end if;                 end do: ## Calculate integrals for RHS vector                for i to m do                F[i] := evalf(Int(theta_init*Efun[i], z = 0 .. d));                if not F[i]::complex(numeric) then                   F[i] := evalf[15](Int(theta_init*Efun[i], z = 0 .. d,                                      digits=15, method=_Dexp, epsilon=1e-12));                   if not F[i]::complex(numeric) then                      print("C",[theta_init*Efun[i], z = 0 .. d]);                   end if;                end if;                end do: ## Create matrix A and RHS vector B                A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):                B := Vector([seq(F[i], i = 1 .. m)]): ## Calculate solve A*x=B               r := LinearSolve(A,B); ## Define Theta(z,t)              theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m): ## Compute v_n for n times constant              for i to m do              v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):              end do: ## Compute v(z,t) from initial conditions              for i to m do              Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):              end do: ## Define v(z,t)              v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m): ##              minp:=min( seq( Re(p[i]), i=1..m) ):              member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'): ## Return all the plots                  return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, H, H__a;                  end proc:
 > # Run the calculation for supplied value of 'xi'
 > # Put the returned quantities in a simple list
 > ans:=CodeTools:-Usage([doCalc(-0.06, 0.001)]): evalf(ans[9]);  ## H value evalf(ans[10]);  ## H bar
 memory used=451.31MiB, alloc change=336.20MiB, cpu time=2.84s, real time=3.06s, gc time=250.00ms
 (1)
 > with(plots):
 > d:= 0.2e-3:
 > testproc := proc(j, u, k) option remember;    local calcvals;    if not [j,u,k]::[numeric,numeric,posint] then       return 'procname'(args);    end if;    calcvals:=doCalc(j,u);    evalf(calcvals[k]); end proc:
 > NN := 8;
 (2)
 > (gridji,rngj,rngi):=[51,101],0..3,-2.0..-0.0; PP[gridji,4]:=plot3d(max(-1e2,testproc(i,j,4)), j=rngj, i=rngi, grid=gridji): PP[gridji,7]:= plot3d(Im(testproc(i,j,7)), j=rngj, i=rngi, grid=gridji):
 (3)
 > ## Surface plot using 3d PP[gridji,4];
 > ## This is the pmin values data[gridji,4]:=op([1,3],PP[gridji,4]):
 > ## Convert Array to Matrix to allow the use listcontplot ConMatrix := Matrix(data[gridji,4]);
 (4)
 > ## Assign the min and max value of pmin to minv and maxv #(minv,maxv):=[min,max](data[gridji,4])[];
 > ## GRB color and the contours (color1,color2):="Red","Yellow": conts:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv];
 >
 (5)
 > colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),                                                i=1..N)])([ColorTools:-Color(color2)[]],                                                          [ColorTools:-Color(color1)[]],                                                          nops(conts))):
 > ##
 > ## Use the Conmatrix and the conts to plot contours and transform the axes subsindets(plots:-listcontplot(ConMatrix, contours=conts,filledregions,thickness=0),            specfunc(anything,THICKNESS),u->THICKNESS(0.2)): LCP[gridji,4]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),                                             lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%): display(LCP[gridji,4],         seq(plot(-2.0,1..1,legend=evalf[4](conts[i]),thickness=15,                  color=colorlist[nops(conts)-i+1]),             i=1..nops(conts)),         size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);
 > ## This deals with the Imaginary part
 > ## Surface plot using plot3d for Imginary part PP[gridji,7]:
 > ##
 > data[gridji,7]:=op([1,3],PP[gridji,7]):
 > ## Convert Array to matrix AMatrix := Matrix(data[gridji,7]);
 (6)
 > (minv,maxv):=[min,max](data[gridji,7])[];
 (7)
 > (color1,color2):="Red","Yellow": conts1:=[minv, seq(minv+(i-1)*(maxv-minv)/(NN+2-1),i=2..NN+1), maxv]/~1.5;
 (8)
 > colorlist:=ColorTools:-Color~(((c1,c2,N)->[seq(c1+(i-1)*(c2-c1)/(N-1),                                                i=1..N)])([ColorTools:-Color(color2)[]],                                                          [ColorTools:-Color(color1)[]],                                                          nops(conts1))):
 > subsindets(plots:-listcontplot(AMatrix,                                contours=conts1-~1e-9,filledregions,thickness=0),            specfunc(anything,THICKNESS),u->THICKNESS(0.2)): LCP[gridji,7]:=plottools:-transform((x,y)->[lhs(rngj)+(rhs(rngj)-lhs(rngj))*(x-1)/(gridji[1]-1),                                             lhs(rngi)+(rhs(rngi)-lhs(rngi))*(y-1)/(gridji[2]-1)])(%): display(LCP[gridji,7],         seq(plot(-2.0,1..1,legend=evalf[4](conts1[i]),thickness=15,                  color=colorlist[nops(conts)-i+1]),             i=1..nops(conts1)),         size=[500,400],legendstyle=[location=right], axes=boxed, labels = ["u", "xi"], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);
 >

 >

@dharr  Thank you...

Thank you

Thank you very much Tomleslie

@Carl Love

Thank you.

Thank you very much for this.

Oh, I see, great

Why you use subsindets??...

@acer Please why do use this line of code?

PC:=subsindets(%,specfunc(anything,THICKNESS),u->THICKNESS(0.2)):

Is the coarse grid because of the -1e2 in the code below?

CodeTools:- Usage(contourplot(max(-1e2,testproc(u,xi,4)), xi=rngxi, u=rngu, grid=griduxi,
contours=conts, thickness=0, coloring=[color1,color2],
filledregions, axes=boxed)):

Appreciation...

@acer Thank you very much for your assistance.

 1 2 3 4 5 6 7 Page 3 of 7
﻿