Dear Maple Primes,

could you, please, help me with numeric integration? I’m new in numeric integration and can’t reach desired precision of a result.
Here is the integral f(x_max) that I try to compute for different values of x_max from the interval 0.025..0.24 :

f:=(xmax)->Int(K*F*Int(G*F,x=x..xmax,method=integrationmethod),x=x0..xmax,method=integrationmethod)

where x₀ is lower limit of outer integral, x₀ := 0.025

and K, F and G are functions of x

K:=x-x0

F:=(a₁+a₂*x+a₃*x²+a₄*x⁵)/(b₁*x+b₂*x²+b₃*x⁶)

G:=exp(c₁+c₂*x+c₃*x⁷)

with

a₁:=8e3; a₂:=6e4; a₃:=3e4; a₄:=1.8e8;
b₁:=9.2e17; b₂:=1.1e18; b₃:=4.6e21;
c₁:=8.202046; c₂:=-12.31377; c₃:=-818043.42;

Please, notice, that G (as well as G*F) is a steeply decreasing function on the interval x = 0.025..0.24.

I get "a seemingly correct" result (that means that f increases as x_max intreases), when I try to plot f(x_max) for the following "guessed" options

Digits:=15
integrationmethod:=_d01akc
plot(f,0.21..0.24,color=black)

What is puzzling me is that I get a different "seemingly correct" result, when I modify the integral f by,
at fist, multiplying G by a constant (for example Const:=1e20; G:=Const*exp(c₁+c₂*x+c₃*x⁷) )
and, second, plotting the f divided by this constant:

plot(f/Const,0.21..0.24,color=red)

The following Figure presents the values of f plotted versus x_max with (red curve) and without (black curve) using of the constant Const:

Dear Primes, could you, please, comment on this difference? Because the only indicator that I have (from the analysis of G, F and K) is that f must be a monotonically (and stricktly) increasing function of x_max.

Please, find the maple worksheet in attachment.

Thank you in advance!
Maks

for_primes_numeric_integration_v02.mw

this is a terrible maple crash just for a simple option !!!

only when i want to print the "x^2+y^2" in the caption of plot it raise exception !
but for some formulae there is no problem !
some technical guide or advise is needed ?
may someone check maple 2018 for this issue too ?

shekofte003.mw

    Hello, colleagues!

    I want to find zero function given by the process that integrates the ODE, but get this warning and do not know what to do with it.  I will be grateful to any advice!!!

T0 := 300; G1 := 1.2*10^6; G2 := 1.2*10^6; k := 1.79333595*10^(-6); l := 1/(12*10^4); U_max := 5.0; th0 := .5; q0 := 0.1e-2:

u(x):=piecewise(x<=-1,-U_max,x>1,U_max,x<=1 and x>-1,0): 

sys := diff(x1(t), t) = x2(t),
           diff(x2(t), t) = l*u(x4(t))-k*sin(2*x1(t)),

           diff(x4(t), t) = -x3(t),

            diff(x3(t), t) = 2*k*cos(2*x1(t))*x4(t);


init:=x1(0)=th0,   x2(0)=q0: 

h1 := proc (t1, alpha_1, alpha_2, th_0, q_0, g1)::float;
local X1, X3, Res;
ptions operator, arrow;
global sys, u;
X1 := eval(x1(t), dsolve({sys, x1(0) = th_0, x2(0) = q_0, x3(0) = alpha_1, x4(0) = alpha_2}, numeric, range = 0 .. t1, output = listprocedure));
X3 := eval(x3(t), dsolve({sys, x1(0) = th_0, x2(0) = q_0, x3(0) = alpha_1, x4(0) = alpha_2}, numeric, range = 0 .. t1, output = listprocedure));
es := X3(t1)+2*g1*X1(t1);
return Res end proc

h2 := proc (t1, alpha_1, alpha_2, th_0, q_0, g2)::float;
local X2, X4, Res;
options operator, arrow;
global sys, u;
X2 := eval(x2(t), dsolve({sys, x1(0) = th_0, x2(0) = q_0, x3(0) = alpha_1, x4(0) = alpha_2}, numeric, range = 0 .. t1, output = listprocedure));
X4 := eval(x4(t), dsolve({sys, x1(0) = th_0, x2(0) = q_0, x3(0) = alpha_1, x4(0) = alpha_2}, numeric, range = 0 .. t1, output = listprocedure));
Res := X4(t1)+2*g2*X2(t1);
return Res
end proc

Test_h1 := proc (x, y)::float;
h1(10.0, x, y, 0.209327944398961e-2, -0.417641536045032e-3, G2)
end proc;

Test_h2 := proc (x, y)::float;
h2(10.0, x, y, 0.209327944398961e-2, -0.417641536045032e-3, G2)
end proc

fsolve([Test_h1, Test_h2]);
Warning, unable to store _EnvDSNumericSaveDigits when datatype=integer[8]
Test_2.mw

I am trying to show visually how many Lie derivatives of two different objects are needed to get a unique solution to a problem, so i want to create a graph of the form:


for the elements of this workseet:
3d_plot_of_Lie_derivatives_against_numelems.mw

I using Maple 18 (not Maple 2018) and I'm trying to figure out how to grab earthquake data from earthquakescanada database from here http://www.earthquakescanada.nrcan.gc.ca//stndon/NEDB-BNDS/bull-en.php using the HTTP requests.

First I used the default search within the web browser, and get a new address which I enter as the URL

HTTP:-Get("http://www.earthquakescanada.nrcan.gc.ca//stndon/NEDB-BNDS/bull-en.php?time_start=2018%2F04%2F23+22%3A47%3A00&time_end=2018%2F05%2F23+22%3A47%3A00&depth_min=0&depth_max=100&mag_min=-3&mag_max=9.9&shape_type=region&radius_center_lat=50&radius_center_lon=-95&radius_radius=1000&region_north=90&region_south=41&region_east=-40&region_west=-150&eq_type_L=1&display_list=1&list_sort=date&list_order=a&tpl_output=html&submited=1"

It takes a long time to download the information and would require HTML surgery but changing the option for output to txt or csv, it's faster and in a much more readable form.  However it's not in a table or Array format, it has become a string.

Is there any way to use ImportMatrix, or ImportData to get a better format of the information?  - both give errors in Maple18.  Or am I stuck trying to use string surgery in Maple 18?  The Import command isn't available until Maple 2016 (I don't mean the Import command within ExcelTools) and I believe that works in Maple 2018 however I'm at a loss for trying to use it in Maple 18. 

Just like the title described, I have encountered an error when I use the command "expand". Actually, I just follow the example, but it doesn't work. Please help me or tell me how can I solve it in other commands.

restart;
alias(epsilon = e, omega = w, omega[0] = w0, t[1] = t1, t[2] = t2); e := proc (t1, t2) options operator, arrow; e end proc; w0 := proc (t1, t2) options operator, arrow; w0 end proc; a := proc (t1, t2) options operator, arrow; a end proc; f := proc (t1, t2) options operator, arrow; f end proc; mu := proc (t1, t2) options operator, arrow; mu end proc;
ode := (D@@2)(u)+2*mu*e*D(u)+w0^2*u+e*w0^2*u^3-e*f*cos(omega*t) = 0;
                                               2  
     @@(D, 2)(u) + 2 mu epsilon D(u) + omega[0]  u

                          2  3                             
        + epsilon omega[0]  u  - epsilon f cos(omega t) = 0
e_oredr := 1;
ode := simplify(subs(D = sum('e^(i-1)*D[i]', 'i' = 1 .. e_oredr+1), ode), {e^(e_oredr+1) = 0});
 / 3         2                                                 
 \u  omega[0]  + 2 (epsilon D[2] + D[1])(u) mu - cos(omega t) f

                  \                   2                   
    + 2 D[1, 2](u)/ epsilon + omega[0]  u + D[1, 1](u) = 0
simplify(collect(%, e), {e^(e_oredr+1) = 0});

u := sum('v[i]*e^i', 'i' = 0 .. e_oredr);
                      epsilon v[1] + v[0]
ode := simplify(collect(ode, e), {e^2 = 0});
for i from 0 to e_oredr do eq[i] := coeff(lhs(ode), e, i) = 0 end do;
                       2                         
               omega[0]  v[0] + D[1, 1](v[0]) = 0
       3         2           2                       
   v[0]  omega[0]  + omega[0]  v[1] + 2 D[1](v[0]) mu

      - cos(omega t) f + 2 D[1, 2](v[0]) + D[1, 1](v[1]) = 0
remove(has, lhs(eq[1]), cos); convert(%(t1, t2), diff);
eq[1] := %-convert(f*cos(sigma*t2+t1*w0), 'exp');

v[0] := A(t2)*cos(w0*t1+B(t2)); convert(%, 'exp'); v[0] := unapply(%, t1, t2);
                         /1                             
       (t1, t2) -> A(t2) |- exp(I (omega[0] t1 + B(t2)))
                         \2                             

            1                              \
          + - exp(-I (omega[0] t1 + B(t2)))|
            2                              /

expand(eq[1]);
Error, (in property/ConvertProperty) invalid input: PropRange uses a 2nd argument, b, which is missing
collect(%, exp(I*w0*t1));
Error, (in collect) invalid 1st argument proc (t1, t2) options operator, arrow; A(t2)*((1/2)*exp(I*(w0*t1+B(t2)))+(1/2)*exp(-I*(w0*t1+B(t2)))) end proc
coeff(%, exp(I*w0*t1));
map(proc (x) options operator, arrow; x*exp(-I*B(t2)) end proc, %);
combine(%, 'exp');
subs(I*B(t2) = I*sigma*t2-I*C(t2), B(t2) = sigma*t2-C(t2), %);
conds := combine(%, 'exp');
                               0

how to change the ode polt color for the below code

with(plots); 

R1 := .1; R0 := .1; m := .1; a := .1; Ha := .1; Nt := .1; Nb := .1; Pr := 6.2; Le := .6; Bi := 1; Ec := .1; k := 1; r := .1; A := 1; 

fcns := {C(y), T(y), U(y), W(y)};

 sys := diff(U(y), `$`(y, 2))+(R1*(diff(U(y), y))-2*R0*W(y))*exp(a*T(y))-a*(diff(U(y), y))*(diff(T(y), y))-Ha = 0, diff(W(y), `$`(y, 2))+(R1*(diff(W(y), y))+2*R0*U(y))*exp(a*T(y))-a*(diff(W(y), y))*(diff(T(y), y))-Ha = 0, diff(T(y), `$`(y, 2))+R1*Pr*(diff(T(y), y))+Pr*Ec*exp(-a*T(y))*((diff(U(y), y))*(diff(U(y), y))+(diff(W(y), y))*(diff(W(y), y)))+Pr*Ha*Ec = 0, diff(C(y), `$`(y, 2))+Pr*Le*R1*(diff(C(y), y))+Nt*(diff(C(y), `$`(y, 2)))/Nb = 0;

bc := U(0) = 0, W(0) = 0, C(0) = 0, (D(T))(0) = Bi*(T(0)-1), U(1) = 0, W(1) = 0, C(1) = 1, T(1) = 0; 

L := [.5, 1.0, 1.5, 2.0]; AP := NULL; 

for k to 4 do R := dsolve(eval({bc, sys}, Ha = L[k]), fcns, type = numeric, AP); AP := approxsoln = R; p1u[k] := odeplot(R, [y, U(y)], 0 .. 1, numpoints = 100, labels = ["y", "U"], style = line, color = ["black", "blue", "red", "pink"]) end do; 

display({p1u[1], p1u[2], p1u[3], p1u[4]})

how to find skin friction value below code

restart

PDEtools[declare]((U, W, T, C)(y), prime = y):

R1 := .1; R0 := .1; m := .1; a := .1; Ha := .1; Nt := .1; Nb := .1; Pr := 6.2; Le := .6; Bi := 1; Ec := .1; k := 1; r := .1; A := 1;

sys := diff(U(y), `$`(y, 2))+(R1*(diff(U(y), y))-2*R0*W(y))*exp(a*T(y))-a*(diff(U(y), y))*(diff(T(y), y))-Ha = 0, diff(W(y), `$`(y, 2))+(R1*(diff(W(y), y))+2*R0*U(y))*exp(a*T(y))-a*(diff(W(y), y)) = 0, diff(T(y), `$`(y, 2))+R1*Pr*(diff(T(y), y))+Pr*Ec*exp(-a*T(y))*((diff(U(y), y))*(diff(U(y), y))+(diff(W(y), y))*(diff(W(y), y)))+Pr*Ha*Ec*((U(y)+m*W(y))*(U(y)+m*W(y))+(W(y)-m*U(y))*(W(y)-m*U(y)))/(m^2+1)^2+Nb*(diff(T(y), y))*(diff(C(y), y))+Nt*(diff(T(y), y))*(diff(T(y), y)) = 0, diff(C(y), `$`(y, 2))+Pr*Le*R1*(diff(C(y), y))+Nt*(diff(C(y), `$`(y, 2)))/Nb = 0:

ba := {sys, C(0) = 0, C(1) = 1, T(1) = 0, U(0) = 0, U(1) = 0, W(0) = 0, W(1) = 0, (D(T))(0) = Bi*(T(0)-1)}:

r1 := dsolve(ba, numeric, output = Array([0., 0.5e-1, .10, .15, .20, .25, .30, .35, .40, .45, .50, .55, .60, .65, .70, .75, .80, .85, .90, .95, 1.00])):

with(plots); 

p1u := odeplot(r1, [y, U(y)], 0 .. 1, numpoints = 100, labels = ["y", "U"], style = line, color = green); 

plots[display]({p1u})

Hi all,

I have the following expression,

where K and K_P are the controller gains, k is the stiffness, b is the damper, m1 and m2 are masses.

How can I find the condition on variables such that the numerator of this expression is greater that zero?

The conditions should appear as inequalities.

Download Routh-Hurwitz_stability.mw

I have four spindle tori governed by the following equations:

I could plot the surfaces using implicitplot3d and I can imagine the volume common to these surfaces but I could not visualize it. So, am looking for a way to plot the volume covered by the surfaces such that f1<0, f2<0, f3<0 and f4<0. I know that it's easy in case of a 2D filled plot but is there any way this could be done for the 3D case? Any mathematical advice as to how to characterize or calculate this volume would also be great.

Using Maple's native syntax, we can calculate the components of acceleration. That is, the tangent and normal scalar component with its respective units of measure. Now the difficult calculations were in the past because with Maple we solved it and we concentrated on the interpretation of the results for engineering. In spanish.

Calculo_Componentes_Aceleracion_Curvilínea.mw

Uso_de_comandos_y_operadores_para_calculos_de_componentes_de_la_aceleración.mw

Lenin Araujo Castillo

Ambassador of Maple

I want to find the i-th root of this equation s := BesselJ(1, alpha)+a1*BesselY(1, alpha)+b1*BeselJ(1, alpha)+c1*BesselK(1, alpha)

When I use the Nextzero it doesnt work

Thanks

Hi,

I couldn't find a command or tool for searching a value in a multi-dimensional list/array/matrix. Does it exist?. 

I wrote this:

where

pSearch: proc for searching

A: list, array, vector or matrix

Examples with a list:

> ss := [[1, 2, 3], [2, 3, 5], [2, 4, 5]]
> j:=1:
># search first index i such that ss[i, j]=2

> nPos := 

# search next index i such that ss[i, j]=2

> nPos := 

But I am not sure if this procedure is efficient for large lists or arrays and I would suggest (if this tool does not exist yet) a more efficient one.

César Lozada

Pretty much the title.

At first I thought something's wrong with my vectors, but even something like <1,1,1>.<2,2,2> produces this error.

I have restart in the very beginning of my worksheet and with(VectorCalculus): right after.

I'm using the 64 bit version of Maple 18 on Windows 10, if this matters.

how can save of resultes as matrix in file.dat?