C_R

3762 Reputation

21 Badges

6 years, 235 days

MaplePrimes Activity


These are answers submitted by C_R

Technically, the injection at t=1 is called an event. An event normally requires reconfiuration of an ODE system. More on dsolve events under ?dsolve,numeric,Events. My understanding of the event options (maybe I am wrong) is that your event (ifelse(t=1,Dose,0)) cannot be given as an option. What you could do instead is splitting the problem at the event into two problems (one before and one after the event) and then plot everything with plots,display. By doing so the event becomes an initial condition for the second problem

However, before you do so, have a look at the event. As you defined it, it does not add any dose to the system because the timespan is infinite small. You need to define an injection amount. I have done this below by introducing a time span epsilon

restart;
k12 := 0.0452821;
k21 := 0.0641682;
k1x := 0.00426118;
Dose := 484;
Inj := Heaviside(t - 1) - Heaviside(t - 1 - epsilon);
epsilon := 0.1;
ode_sys := diff(SA2(t), t) = SA1(t)*k12 - SA2(t)*k21, diff(SA1(t), t) = Inj - k1x*SA1(t);
ics := SA2(0) = 0, SA1(0) = 0;
Solution1 := dsolve({ics, ode_sys}, {SA1(t), SA2(t)}, numeric);
plots:-odeplot(Solution1, [[t, SA1(t), color = "red", thickness = 1], [t, SA2(t), color = "blue", thickness = 1]], t = 0 .. 500, labels = ["Time t", "Salicylic acid"], labeldirections = [horizontal, vertical]);

 

Why not remove the differential from the start?

In your example you could have devided expression (1) by dz. This effectively transforms (1) into a differential quotient.

Maple does not accept differentials for integration but can work this way with differential quotients.

restart

NULL

`τ_` = int(diff(`τ_`(z), z), z = l .. -l)

`τ_` = int(diff(`τ_`(z), z), z = l .. -l)

(1)

diff(`τ_`(z), z) = -I*z*B__0*x^2*_j/d^2-I*y*B__0*x^2*_k/d^2

diff(`τ_`(z), z) = -I*z*B__0*x^2*_j/d^2-I*y*B__0*x^2*_k/d^2

(2)

subs(diff(`τ_`(z), z) = -I*z*B__0*x^2*_j/d^2-I*y*B__0*x^2*_k/d^2, `τ_` = int(diff(`τ_`(z), z), z = l .. -l))

`τ_` = int(-I*z*B__0*x^2*_j/d^2-I*y*B__0*x^2*_k/d^2, z = l .. -l)

(3)

value(%)

`τ_` = (2*I)*y*B__0*x^2*_k*l/d^2

(4)

NULL


 Mathematically, this is the exact way and (I assume) the reason why we don't find differentials in Maple, even though they are still widely used in physics and engineering.

Download RemoveDifferential_wothout_differential.mw

with 2024

(normal@expand@combine)((2), trig)

same as

normal(expand(combine(expr2, trig)))

 

You can use format strings to format the output as you want

a := 3*89/100.;
b := 3*21/100.;
printf("%g,%g,%6.2f,%9.2f", a, b, a, b);
                        a := 2.670000000

                       b := 0.6300000000

2.67,0.63,  2.67,     0.63

Newer versions of Maple have the context pannel (or right click on the output) for numeric formating. I don't know if that is possible in Maple 2016.

 

Prefixnotation is especially useful in functional programming where a sequence of commands is grouped using the composition operator @. It is way of separating functions (i.e. Maple commands) from arguments. Otherwise, nested commands with allot of parentheses have to be used.

There is little documentation on that and it is IMO not self-explaining. I got a pretty good explanation here in MaplePrimes. And here is a nice example of a self-made “unary prefix operator” that simplifies a rational expression. (I have put operator in quotation marks since depending on what is does it can equally be called function or command in Maple. In this example it is rather a command.)

I agree that for only two operands infix notation is often preferable. But in functional programming the number of operands (or better arguments if interpreted from a procedure perspective) can change in the course of processing composed commands.

I think a bit more explanatory documentation on top of dharr's examples would not harm.

I hope this completes the picture why Maple sometimes provides infix and prefix notation for the same operation. Functional programming without some prefix variants would be less powerful.

Units_and_Integrals_reply.mw

I have fixed a few errors

I assume that you are familar with Maple styles and in your question you want, for example, to change only the font of all Maple input in an existing document.

What you have to do is:

  • hide all other screen content using view-> show/hide contents
  • select all remaining content with the mouse
  • change the font
  • show the hinden content using view-> show/hide contents

Units in (4) are correct if t in (4) is evaluated with units. With units standard int must be mapped over a column vector (I don't know why that is).

restart

with(VectorCalculus); SetCoordinates('cartesian[x, y, z]')

with(ScientificConstants)

"with(Units[Simple]):"

NULL

g := evalf(Constant(g, units))

9.80665*Units:-Unit(m/s^2)

(1)

`#mover(mi("\`v__0\`"),mo("&rarr;"))` := `<,>`(v[0]*Unit('m'/'s'), 0*Unit('m'/'s'), 0*Unit('m'/'s'))

Vector[column](%id = 36893490633426056244)

(2)

`#mover(mi("a"),mo("&rarr;"))` := `<,>`(0*Unit('m'/'s'^2), g, 0*Unit('m'/'s'^2))

Vector[column](%id = 36893490633436725244)

(3)

int(`#mover(mi("a"),mo("&rarr;"))`, t)

Vector[column](%id = 36893490633436713196)

(4)

"eval(?,t=2&lobrk;s&robrk;)"

Vector[column](%id = 36893490633436699932)

(5)

restart;
with(Units:-Standard):

`#mover(mi("a"),mo("&rarr;"))` := <0*Unit(('m')/'s'^2), g, 0*Unit(('m')/'s'^2)>

Vector(3, {(1) = 0, (2) = g, (3) = 0})

(6)

map(int,(6),t)

Vector(3, {(1) = 0, (2) = g*t, (3) = 0})

(7)

int~((6),t)

Vector(3, {(1) = 0, (2) = g*t, (3) = 0})

(8)

Download Units_and_Integrals_reply.mw

I assume that you are able to generate plots for the exact and the approximate solution using the plot command.

Each of these plots have to be assigned to a name with the assignment operator :=

Then you can plot boths solutions in one plot using plots:- display command.

If you have these solutions in a worksheet and you can't manage to compare them, you could upload them using the green arrow.

To your question: How to calculate this integral?

My preliminary answer is: Don't use Maples sophisticated algorithms.

Zoomed in it looks like this

This could explain why Maples sophisticated integration methods have difficulties. They will try to follow all the oscillations.

Annother reason could be: The function you want to integrate is not always real valued also (see the attachment).

Could you tell a bit more about your integrand?

PS:: I have just seen Roubens answer. Zoomed it looks also jagged.

Integral_reply.mw

Do you want the new constant to have the highest index? If filling of gaps is acceptable

n:=1; 
do n++ until not(has(_C||n,myconstants)):
new_constant := _C ||(n)
 

Otherwise, in your example myconstants is a set and the element with the highest index is to the right.

n:=1; 
do n++ until is(_C||n=myconstants[-1]):
new_constant := _C ||(n)

I guess there are shorter ways

but you will get three bulky solutions

For clarity I assume in the following that you combine the coefficients into new ones.
Your equations will look like this

eq1:=cos(3*phi)-a*cos(phi)=b

cos(3*phi)-a*cos(phi) = b

(1)

eq2:=sin(3*phi)-a*sin(phi)=c

sin(3*phi)-a*sin(phi) = c

(2)

 

Solve for phi

expand(eq1,trig);
sols:=solve(%,[phi])

4*cos(phi)^3-3*cos(phi)-a*cos(phi) = b

 

[[phi = arccos((1/6)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-6*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3))], [phi = Pi-arccos((1/12)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-3*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+6*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)))], [phi = arccos(-(1/12)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+3*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+6*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)))]]

(3)

You get 3 solutions for for phi that you can substitute one after each other into eq2

NULL

simplify(subs(sols[1],eq2))

(1/18)*(-9*(((2/3)*a-2)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3)+(3*b+(1/3)*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+(a+3)^2)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3))^(1/2)*((-a+3)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3)+(9*b+(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+3*(a+3)^2)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3) = c

(4)

 

Depending on the initial parameters you might be able to simplify further and select an equation that describes your problem

Download elliminate_phi.mw

Maple is capable to integrate up to the first 200 roots without any special methods.
 

evalf(Int(abs(sin(x^4))/(sqrt(x) + x^2), x = 0 .. 5))

 

50000 roots are reached at about x=20.

 

evalf(eval((Pi*i)^(1/4), i = 50000))

However, maxintervals=50000 is in any case not sufficient to place between all roots up to infinity at least one interval. But this is somehow  irrelevant, because Maple does not even get to the point where 50000 intervals are reached.

What happens in your case: Maple subdivides the integrations range applying numerical analysis to the integrand and stops at 200 which matches the 200 roots (see above).

Error, (in evalf/int) NE_QUAD_MAX_SUBDIV:
  The maximum number of subdivisions has
  been reached: max_num_subint = 200

Reducing the numerical accuracy to 2 significant digits does prevent Maple from detecting more than 200 roots. The evaluated integrand becomes too small and is rounded to zero. This limits the detected roots.

Increasing to 3 allows you to integrate up to 17

evalf[3](Int(abs(sin(x^4))/(sqrt(x) + x^2), x = 0 .. 16.999))

Increasing to 4, and you can go up to 13

evalf[4](Int(abs(sin(x^4))/(sqrt(x) + x^2), x = 0 .. 13.0001))

I think this challenging integral needs to be attacked from a different angle. Maybe something with limits or an adaptive method that integrates only up to the next 200 roots until a convergence criteria is reached.

Update:
The error messages the OP was observing were contradictory in the sense that in one instance maxintervals was ignored, in the other not. It is still unclear why forcing maxintervals did not work in one case. Anyway, to substantiate my answer, I have added a worksheet which highlights essential parts of Maples integration control. It can be seen how evalf[n] sets integration error tolerances, at which instance the integrand is set to zero, and how Maple iteratively subdivides an integration interval, tries to apply default methods to a new subinterval, switches to specific methods and changes the methods when numerically required.

adaptive_integration.mw

If your region/domain is any other than the coords option of plot offers, then you have to write your own plot routine. What I mean by this is that in (for example) cartesian coordinates only rectangular regions can be mapped.

In the example below I have punched a hole in the domain of the quadrant Re(z)>0, Im(z)>0.

restart;
f:=z^2;# try f:=z to see how the domain is restricted

z_map:=proc(f,re,im) 
  if(sqrt((re-5)^2+(im-3)^2)>=3)then
    eval(f,z=re+I*im);
  else
    NULL;
  end if;
end proc;

p_re:=plots:-display(seq(plot([Re('z_map(f,re,im)'),Im('z_map(f,re,im)'),im=0..10]),re=0..20)):
p_im:=plots:-display(seq(plot([Re('z_map(f,re,im)'),Im('z_map(f,re,im)'),re=0..20]),im=0..10),color=green):
plots:-display(p_re,p_im,scaling=constrained)

The above could be expanded to a routine where you could pass the restrictions as arguments.

Update:

A triangular domain can be found here.

 

X := Matrix([[1, -X3_3/2 - 1/2, 0, -X2_3], [-X3_3/2 - 1/2, -2*X3_4 - 1, X2_3, 0], [0, X2_3, X3_3, X3_4], [-X2_3, 0, X3_4, 1]]);
vars := [X3_4, X3_3, X2_3];
w := A^3 - A;
rootz := RootOf(w, A);
Pols := [(-A^2 + 1)/(3*A^2 - 1), (-A^2 - 1)/(3*A^2 - 1), A*(3*A^2 - 1)*1/(3*A^2 - 1)];
vals := {allvalues(eval(Pols, A = rootz))};
ecs := map(x -> convert(vars = x, listofequations)[], vals);
Xs := [seq(eval(X, ecs[x]), x = 1 .. numelems(vals))];
Xz := [seq(subs(ecs[x], X), x = 1 .. numelems(vals))];
with(LinearAlgebra);
Eigs := map(x -> Eigenvalues(x), Xs);
BolEigs := map~(type, Eigs, poszero);
evalf(Eigs);

For BolEigs there is probably a shorter version without seq.

Edit: Must be poszero not positive. I have corrected that in the above

Edit2: BolEigs replaced by something shorter in the above

 

First 7 8 9 10 11 12 13 Last Page 9 of 18