Advanced Search

MaplePrimes Questions

fixed point method...

Asked by: Zeineb 530
rootfinding fsolve
April 25 2022
1  2

Dear all

I have a nonlinear system of algebraic equations, I would like to solve it  without using fsolve and Newton's method. 
Maybe one can use, fixed point method or Broyden's method

Fixed_Broyden_method.mw

Can we solve the system of 4 equations, using the proposed methods

thank you 

using printf to print a palletes using as an outpu...

Asked by: qtngmin 15
printf 2d-output
April 25 2022
2  2

Hi everyone, I have some questions about the printf func in maple.

As an image I attach below, I want to print this to the console so that I can presentation it for my teacher, because I'm using this palletes many times, only showing the result will make me and my teacher hard to check the right result. I trying using the sum palletes as a tring ("%s") and also the result as ("%f") too, but it wont work. Please help me with this, thanks a lot

Solve & Trig: Why is the output different when arc...

Asked by: C_R 2592
solve arctan inverse-kinematics
April 25 2022
2  8

Solve produces different output in the attachment depending on how it is used. Why is that and how can simplification to arctan(y/z) be avoided? Arctan(y/z) only gives correct angles for positive y and z.  I prefer arctan(y,z) output that I can subsequently simplify to the y and z ranges of interest (if possible). Imagine “wrong” simplification of complex algebraic output (e.g., from inverse kinematics).

Arctan.mw

display problem of constant of integration when us...

Asked by: nm 9889
alias latex
April 25 2022
3  1

I am using the method of alias(seq(c[k] = _C||k, k = 1..10)); for better latex of constants generated from solving an ode as recommended. See this for example. 

This works well 99.99% of the time. But now I noticed this in Maple 2022. Is this a display issue? When the constant is inside an inert Int it does not display the same as the other constant outside. Also the Latex is not the same. Even though lprint shows they are both correct.

This is the worksheet itself


 

> 

restart;
> 

interface(version);

`Standard Worksheet Interface, Maple 2022.0, Windows 10, March 8 2022 Build ID 1599809`

> 

alias(seq(c[k] = _C||k, k = 1..10));
ode:=x*diff(y(x),x$2)-cos(x)*diff(y(x),x)+sin(x)*y(x)=2;
sol:=dsolve(ode);
lprint(sol)

c[1], c[2], c[3], c[4], c[5], c[6], c[7], c[8], c[9], c[10]

x*(diff(diff(y(x), x), x))-cos(x)*(diff(y(x), x))+sin(x)*y(x) = 2

y(x) = (c[2]+Int((c[1]+2*x)/(exp(Ci(x))*x^2), x))*exp(Ci(x))*x

y(x) = (c[2]+Int((c[1]+2*x)/exp(Ci(x))/x^2,x))*exp(Ci(x))*x

> 

latex(sol)

y \! \left(x \right) =
\left(c_{2}+\textcolor{gray}{\int}\frac{\mathit{c[1]} +2 x}{{\mathrm e}^{\mathrm{Ci}\left(x \right)} x^{2}}\textcolor{gray}{d}x \right) {\mathrm e}^{\mathrm{Ci}\left(x \right)} x

> 

restart;
> 

ode:=x*diff(y(x),x$2)-cos(x)*diff(y(x),x)+sin(x)*y(x)=2;
sol:=dsolve(ode);
lprint(sol)

x*(diff(diff(y(x), x), x))-cos(x)*(diff(y(x), x))+sin(x)*y(x) = 2

y(x) = (_C2+Int((_C1+2*x)/(exp(Ci(x))*x^2), x))*exp(Ci(x))*x

y(x) = (_C2+Int((_C1+2*x)/exp(Ci(x))/x^2,x))*exp(Ci(x))*x

> 

 


 

Download april_25_2022.mw

 

 

animate plot in Maple...

Asked by: charlie_fcl 110
animation plotting
April 24 2022
1  5

Hi,

I'm experimenting with some plot functions in Maple and cannot make the animate function to work, even for the help example (attached). It works fine when I insert it here but not in my worksheet.

with(plots)

NULL

plots[animate](plot, [A*sin(x), x = 0 .. 10], A = 0 .. 2)

  

NULL

Download animate_prob.mw

Any ideas?

loops in MapleFlow...

Asked by: charlie_fcl 110
do-loop
April 24 2022
3  2

Hi,

I like to know if it is possible to use a for loops in MapleFlow. I've tried a couple of times but it didn't work. For example,

n:=6
a:=Matrix(n,n)

for l from 1 to n do

   for c from 1 to n do

     a[l,c]:=2^(l+c)

   end do

end do

a=

 

doesn't produce any result for a

Best regards,

Numerical Scheme for Integro-differential equation...

Asked by: Tamour_Zubair 80
numeric differential-equations integro-differential-equation
April 23 2022
1  1

Hello;

Hope you are fine. Can i apply numerical scheme on maple for the following problem. This in integro-differential equation i think. Waiting for kind response.

Thanks

 

"Help!for this series."...

Asked by: Saha 75
solve homework
April 23 2022
0  1


could you please help me ,the maple code for this given series.

> 

restart
> 

U[0](x) := x;

x

 (1)
> 

"U[k+1](x):=solve((k+1)*U[k+1](x)+(x*(-1)^((k-1)/(2)))/(k!)-x^(2)*((e)^(x))/(10){6/(k!)-sum((2^(k[1]))/(k[1!])(5*delta[k[]-k[1]](x)+(2^(k-k[1])*(-1)^((k-k[1])/(2)))/((k-k[1])!)+(2^(k-k[1]+1)*(-1)^((k-k[1]-1)/(2)))/((k-k[1])!)),k[1]=0..k)}-(cos^(2)(x)+sin^()(x))*((∂)^2)/((∂)^( )x^2) [U[k](x)]-(e)^(x)[sum(1/(k[1]!){1/(k-k[1])(sum(sum(1/(k[3]!)*U[k[2]-k[3]](x)*U[k-k[1]-k[2]-1](x)},k[1]=0..k),k[2]=0..k-k[1]-1),k[3]=0..k[2])),U[k+1](x)]; od;"

Error, unable to match delimiters

"U[k+1](x):=solve((k+1)*U[k+1](x)+(x*(-1)^((k-1)/2))/(k!)-x^2*((e)^x)/10{6/(k!)-sum((2^(k[1]))/(k[1!])(5*delta[k-k[1]](x)+(2^(k-k[1])*(-1)^((k-k[1])/2))/((k-k[1])!)+(2^(k-k[1]+1)*(-1)^((k-k[1]-1)/2))/((k-k[1])!)),k[1]=0..k)}-(cos^2(x)+sin(x))*((∂)^2)/(∂x^2) [U[k](x)]-(e)^x[sum(1/(k[1]!){1/(k-k[1])(sum(sum(1/(k[3]!)*U[k[2]-k[3]](x)*U[k-k[1]-k[2]-1](x)},k[1]=0..k),k[2]=0..k-k[1]-1),k[3]=0..k[2])),U[k+1](x)]; od;"

  
> 

``


 

Download Chapter_6-Example-6.5.4.mw

How to make animation of matrix plotting ?...

Asked by: PavelGR 5
animation plotting incomplete-question
April 23 2022
1  2

How to make the matrix construction by points animated?
I sketched this code - it doesn't work(((
 

a := Matrix([stoimostyM, izmeniya]);
plots[animate](plot, [a[n, m]], n = 1 .. 2970, m = 1 .. 2, numpoints = 50, frames = 100);
Error, bad index into Matrix

Belousov–Zhabotinsky equation ...

Asked by: Pepini 255
animation differential-equations laplace
April 23 2022
2  2

Hi! Do you know maybe how to solve equation with Laplace operator in Maple like BZ equation?

restart;
a := 0.75;
rho := u(t) + v(t) + w(t);
ode := diff(u(t), t) = 10*Delta(u(t)) + u*(-a*v - rho + 1);
ode1 := diff(v(t), t) = 10*Delta(v(t)) + v*(-a*w - rho + 1);
ode2 := diff(w(t), t) = 10*Delta(w(t)) + w*(-a*u - rho + 1);

 

Edit: Sorry I guess it should be function of three variables so u,v,w depends on (x,y,z) not strictle from time

I am wondering how to animate something like this from BZ equation:

How can I generate the following list from L=[x,y,...

Asked by: Chorux 60
combination list
April 23 2022
1  5

Given a list L=[x,y,z], what is the best way to generate the following sequential expression?
L1=[[x],[y],[z],[x,x],[x,y],[[x,z],[y,y],[y,z],[z,z],[x,x,x],[x,x,y],[x,x,z],[x,y,y],[x,y,z],[y,y,y],[y,y,z],[y,z,z],[z,z,z]].

With L1 in mind (a list of list of at most three variables), how can one generate a list of lists of any number of variable from L.

I am currently out of options, a response will be highly appreciated.

Thanks

ERROR: Flow does not evaluate to a module...

Asked by: rs4670 10
April 22 2022
2  1

 

I have Maple 2021 and Mapleflow 2021 installed. The Maple 2021 is working fine. However, MapleFlow gives the error as follows. ERROR: Flow does not evaluate to a module. My system OS is Windows 10 updates till April 2022

 

 

Will anyone be able to speed up this code even mor...

Asked by: Dr. Venkat Subramanian 411
sparse performance linearalgebra
April 22 2022
3  17

Details about our approach for phase-field models are provided in a paper just submitted. 
https://ecsarxiv.org/k2vu6/
 

Maple code (based on UMFPACK linearsolver and different compiled procedures for marching in time) is given here. Run the code for small values of NN and MM and keep increasing them. The goal is to do 100, 200, 400, etc and still get the code to run very fast.
 

> 

restart:
> 

t11:=time():
> 

t12:=time[real]():
> 

Digits:=15;

15

 (1)

NN is the number of node points (elements) in the X and MM is the number of elements in the Y direction. delta is the applied current density at the top (Y =1). tf is the final time for simulation. vel is the velocity constant v in the paper. ki0 is the scaled exchange current density k in the paper. This code can be run for positive values of delta. This simulates plating. At the end of simulation, changing delta to negative values and rerunning the code will automatically used the geometry at the end of plating.
Ydatstore stores the geometry at every point in time. Phiaveadd stores the total liquid phase in the domain at any point in time.
Users can change NN, delta, tf, vel, ki0, MM just in this line and choose Edit execute worksheet to run for different design parameters.

Users can modify the call for y0proc for choosing different models.

Users can modify the call for HF to run first-order upwind, ENO2 or WENO3 methods. NN and MM should be even numbers.

> 

NN:=100;delta:=0.1;tf:=2.0;vel:=1.0;ki0:=1.0;MM:=NN;

100

  

.1

  

2.0

  

1.0

  

1.0

  

100

 (2)
> 

N:=NN+2:h:=1.0/NN:M:=MM+2;Ny:=MM;

102

  

100

 (3)
> 

 

Initial geometry, Model 1, semicircle in a square

> 

y0proc1:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
> 

N:=NN+2:h:=1.0/NN:w:=h/2:
> 

M:=MM+2;Ny:=MM;
> 

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=xx^2+yy^2;
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
> 

 

Model 2, square in a square

> 

y0proc2:=proc(NN,MM,Y00)# square inside a circle, Model 2
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
> 

N:=NN+2:h:=1.0/NN:w:=h/2:
> 

M:=MM+2;Ny:=MM;
> 

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if xx <=0.3 and yy<=0.3 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
> 

 

> 

 

Initial geometry, Model 3, electrodeposition problem trenches and via

> 

y0proc3:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
> 

N:=NN+2:h:=1.0/NN:w:=h/2:
> 

M:=MM+2;Ny:=MM;
> 

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if abs(xx-0.5) >0.2 and yy<=0.5 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
> 

 

> 

 

Initial geometry, Model 4, Gaussian Seed at the bottom

> 

y0proc4:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
> 

N:=NN+2:h:=1.0/NN:w:=h/2:
> 

M:=MM+2;Ny:=MM;
> 

for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=0.1+0.1*exp(-500.*(xx-0.5)^2);
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((yy-rr)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
> 

 

> 

y0proc:=eval(y0proc1):#choose different models using y0proc2, etc.
> 

Y00:=Matrix(1..N,1..M,datatype=float[8]):
> 

evalhf(y0proc(NN,MM,Y00)):
> 

if delta<0 then read("Y0data.m"):end:
> 

if delta<0 then read("tdata.m"):end:
> 

tf;

2.0

 (4)
> 

p0:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):
> 

printlevel:=2;

2

 (5)
> 

eq:=Array(1..N,1..M):
> 

Nx:=NN:
> 

 

> 

 

Next, boundary conditiosn at X  = 0, X =1, Y = 0, Y = 1 are specified below, but these equations are not used and optimally coded inside the procedure Eqs11.

> 

for i from 1 to 1 do for j from 1 to M do eq[i,j]:=-Phi[i,j]+Phi[i+1,j]:od:od:
> 

for i from N to N do for j from 1 to M do eq[i,j]:=Phi[i-1,j]-Phi[i,j]:od:od:
> 

for i from 1 to N do for j from 1 to 1 do eq[i,j]:=Phi[i,j+1]-Phi[i,j]:od:od:
> 

for i from 1 to N do for j from M to M do eq[i,j]:=Phi[i,j-1]-Phi[i,j]+delta*h:od:od:
> 

Eqs:=Vector(N*(M)):
> 

N,M;

102, 102

 (6)
> 

 

Residues at different points in X and Y are coded in the Eqs11 procedure. Y0 is the input phase-field parameter (2D Matrix). The potential y is a vector to expedite the calculation of residue.

> 

Eqs11:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,ff::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):

for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
ff[i1]:=
(Y0[i,j]+Y0[i,j+1])*(y[i1+N]-y[i1])
-(Y0[i,j]+Y0[i,j-1])*(y[i1]-y[i1-N])
+(Y0[i+1,j]+Y0[i,j])*(y[i1+1]-y[i1])
-(Y0[i,j]+Y0[i-1,j])*(y[i1]-y[i1-1])
-ki0*y[i1]*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+N]:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:
> 

 

The Jacobian for the residues is coded int he procedure Jac1.

> 

Jac1:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,j00::Matrix(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):

for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
j00[i1,i1]:=-4*Y0[i,j]-Y0[i,j+1]-Y0[i,j-1]-Y0[i+1,j]-Y0[i-1,j]-ki0*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
j00[i1,i1+1]:=Y0[i+1,j]+Y0[i,j]:j00[i1,i1-1]:=Y0[i-1,j]+Y0[i,j]:
j00[i1,i1+N]:=Y0[i,j+1]+Y0[i,j]:j00[i1,i1-N]:=Y0[i,j-1]+Y0[i,j]:
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1+1]:=1.00:
#ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-1]:=1.00:
#ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
j00[i1,i1]:=-1.0:j00[i1,i1+N]:=1.00:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-N]:=1.00:
#ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:
> 

Eqs11:=Compiler:-Compile(Eqs11):
> 

 

> 

ntot:=N*(M);

10404

 (7)
> 

Ntot:=ntot;

10404

 (8)
> 

 

First-order Upwind method is coded in the procedure UpW.

> 

UpW:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8];
#e1:=1e-15:

h:=1.0/(N-2):

vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:


vx1:=(Y00[i,j]-Y00[i-1,j])/h:
vx2:=(Y00[i+1,j]-Y00[i,j])/h:


vy1:=(Y00[i,j]-Y00[i,j-1])/h:
vy2:=(Y00[i,j+1]-Y00[i,j])/h:

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
F0[i,j]:=nx:
od:od:

end proc:
> 

 

Second-order ENO2 method is coded in the procedure ENO2.

> 

ENO2:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],alpha::float[8],beta::float[8];


h:=1.0/(N-2):
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:

sdx:=(Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j])/h:sdy:=(Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1])/h:
vxb:=0:vxf:=0:vyb:=0:vyf:=0:
if i = 2 then sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j])/h:else sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j])/h:end:
if sdx*sdxb>=0 then s1x:=1.0 else s1x:=0.0:end:
vx1:=(Y00[i,j]-Y00[i-1,j])/h+0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxb)):

if i = N-1 then sdxf:=(Y00[i+1,j]-2*Y00[i+1,j]+Y00[i,j])/h:else sdxf:=(Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j])/h:end:
if sdx*sdxf>=0 then s1x:=1.0 else s1x:=0.0:end:
vx2:=(Y00[i+1,j]-Y00[i,j])/h-0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxf)):


if j = 2 then sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1])/h:else sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2])/h:end:
if sdy*sdyb>=0 then s1y:=1.0 else s1y:=0.0:end:
vy1:=(Y00[i,j]-Y00[i,j-1])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyb)):

if j = M-1 then sdyf:=(Y00[i,j+1]-2*Y00[i,j+1]+Y00[i,j])/h:else sdyf:=(Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j])/h:end:
if sdy*sdyf>=0 then s1y:=1.0 else s1y:=0.0:end:
vy2:=(Y00[i,j+1]-Y00[i,j])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyf)):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):

F0[i,j]:=nx:
od:od:

end proc:
> 

 

Third-order WENO3 method is coded in the procedure WENO3.

> 

WENO3:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8],e1::float[8];
e1:=1e-6:

h:=1.0/(N-2):

vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
phix:=(Y00[i+1,j]-Y00[i-1,j])/2/h:
phiy:=(Y00[i,j+1]-Y00[i,j-1])/2/h:

if i = 2 then sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j]: else sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j]:end:
sd:=Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j]:
if i = N-1 then sdf:=Y00[i,j]-2*Y00[i+1,j]+Y00[i+1,j]: else sdf:=Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):
vx1:=phix-0.5*w1/h*(sd-sdb):
vx2:=phix-0.5*w2/h*(sdf-sd):

if j = 2 then sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1]:else sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2]:end:
sd:=Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1]:
if j = M-1 then sdf:=Y00[i,j]-2*Y00[i,j+1]+Y00[i,j+1]: else sdf:=Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):

vy1:=phiy-0.5*w1/h*(sd-sdb):
vy2:=phiy-0.5*w2/h*(sdf-sd):

if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:

if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
#nx:=sqrt(vx1^2+vx2^2+vy1^2+vy^2):

F0[i,j]:=nx:
od:od:

end proc:
> 

 

> 

UpW:=Compiler:-Compile(UpW):
> 

ENO2:=Compiler:-Compile(ENO2):
> 

WENO3:=Compiler:-Compile(WENO3):
> 

F0:=Matrix(1..N,1..M,datatype=float[8]):
> 

h/0.1/2;

0.500000000000000e-1

 (9)

 

> 

Nx;

100

 (10)
> 

PhiAdd:=proc(N::integer,Phi0::Vector(datatype=float[8]),db::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phi0[i]:=Phi0[i]+db[i]:od:
end proc:
> 

PhiAdd:=Compiler:-Compile(PhiAdd):
> 

Ntot:=ntot:
> 

Phi0:=Vector(1..Ntot,datatype=float[8]):
> 

printlevel:=1:
> 

ff:=copy(Phi0):
> 

 

> 

Phi0:=Vector(1..Ntot,datatype=float[8]):
> 

j00:=Matrix(1..Ntot,1..Ntot,datatype=float[8],storage=sparse):
> 

evalf(Eqs11(Y00,Phi0,delta,ki0,ff));

0.100000000000000e-2

 (11)
> 

Jac1(Y00,Phi0,delta,ki0,j00):
> 

db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
> 

Phi0:=Phi0+db:
> 

V[0]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;

HFloat(0.31581237885872887)

 (12)
> 

TT[0]:=0;tt:=0:

0

 (13)
> 

vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
<