Advanced Search

Maple Questions and Posts

These are Posts and Questions associated with the product, Maple
View unanswered posts
Maple Questions and Posts Feed

How to find common members...

Asked by: brian bovril 914
list intersect
November 06 2024
2  4

Hi. I'm comparing two 7x3-man teams lists. I can see there are at least 2 teams that have common members, new[1] and past[4], new[4] and past[3]. There may be others but the chat gpt code below doesn't definitively find any of them.

common.mw

Library size distribution...

Asked by: Anthrazit 885
library
November 05 2024
3  3

The company library with all the good stuff has passed a 100 Mb size, and I wonder if it is possible to find out, which parts of the library that use most of the storage space?

Nice recursion to solve...

Asked by: Alfred_F 345
sum recurrence limit
November 05 2024
0  2

After exertion with ordinary differential equations now relaxation:

Determine the formation law, limit and sum limit for
u_n+3=(13/12)*u_n+2 - (3/8)*u_n+1 + (1/24)*u_n .
Starting values ​​u_1=0, u_2=1, u_3=1.

Introducing the Student Success Platform

by: Karishma 825 Maple , Maple Calculator , Maple Learn
launch
November 05 2024
6  2

We’re thrilled to announce the launch of our new Student Success Platform! Over the past several months, our academic team has dedicated itself to understanding how we can better support institutions in addressing their concerns around student retention rates. The numbers tell a concerning story: In the U.S., nearly 25% of first-year undergraduates don’t complete their studies, and in STEM fields, the numbers are even higher. In both STEM programs and non-STEM programs with math gateway courses, struggles with math are often a key reason students do not, or cannot, continue their studies. This has a profound impact on both the students’ futures and the institution’s revenue and funding.

From what we’re hearing from institutions and instructors, one of the most pressing issues is the lack of readiness among first-year students, particularly in math courses. With larger class sizes and students arriving with varying levels of preparedness, instructors face challenges in providing the personalized support that is essential. Additionally, many students don’t fully utilize existing resources, such as office hours or TA sessions, which increases their risk of falling behind and ultimately dropping out.

Our new Student Success Platform is designed to tackle these issues head-on. It combines all of our existing tools with exciting new features to help students succeed on their own terms—without adding to instructors' already busy workloads. The early feedback has been fantastic, and we can’t wait for you to see the impact it can make.

You can read more about the Student Success Platform here: https://www.maplesoft.com/student-success-platform/

 

How make a table for numerical results between exa...

Asked by: salim-barzani 1610
table numeric exact
November 05 2024
0  2

If we calculating it take to much time but if we make a procedure it will be more effectable for such example, i want the exact and approximat and error

> 

restart
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(inttrans)
> 

with(SolveTools)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

 (1)
> 

NULL
> 

eq := diff(y(x, t), `$`(t, 2))+(1+x)*(diff(y(x, t), x))-y(x, t) = 2*y(x, t)^3

diff(diff(y(x, t), t), t)+(1+x)*(diff(y(x, t), x))-y(x, t) = 2*y(x, t)^3

 (2)
> 

eq1 := laplace(eq, t, s)

s^2*laplace(y(x, t), t, s)-s*y(x, 0)+laplace(diff(y(x, t), x), t, s)*x-laplace(y(x, t), t, s)-(D[2](y))(x, 0)+laplace(diff(y(x, t), x), t, s) = 2*laplace(y(x, t)^3, t, s)

 (3)
> 

eq2 := subs({y(x, 0) = 1, (D[2](y))(x, 0) = 1}, eq1)

s^2*laplace(y(x, t), t, s)-s+laplace(diff(y(x, t), x), t, s)*x-laplace(y(x, t), t, s)-1+laplace(diff(y(x, t), x), t, s) = 2*laplace(y(x, t)^3, t, s)

 (4)
> 

eq3 := s^2*laplace(y(x, t), t, s) = s-laplace(diff(y(x, t), x), t, s)*x+1+laplace(diff(y(x, t), x), t, s)+2*laplace(y(x, t)^3, t, s)+laplace(y(x, t), t, s)

s^2*laplace(y(x, t), t, s) = s-laplace(diff(y(x, t), x), t, s)*x+1+laplace(diff(y(x, t), x), t, s)+2*laplace(y(x, t)^3, t, s)+laplace(y(x, t), t, s)

 (5)
> 

eq4 := expand(eq3/s^2)

laplace(y(x, t), t, s) = 1/s-laplace(diff(y(x, t), x), t, s)*x/s^2+1/s^2+laplace(diff(y(x, t), x), t, s)/s^2+2*laplace(y(x, t)^3, t, s)/s^2+laplace(y(x, t), t, s)/s^2

 (6)
> 

NULL
> 

"u[0](x):=invlaplace(1/s+1/(s^2),s,x)"

proc (x) options operator, arrow, function_assign; invlaplace(1/s+1/s^2, s, x) end proc

 (7)
> 

u[0](x)

1+x

 (8)
> 

n := N

N

 (9)
> 

k := K

K

 (10)
> 

f := proc (u) options operator, arrow; u^3 end proc

proc (u) options operator, arrow; u^3 end proc

 (11)
> 

for j from 0 to 3 do A[j] := subs(lambda = 0, (diff(f(seq(sum(lambda^i*u[i](x), i = 0 .. 20), m = 1 .. 2)), [`$`(lambda, j)]))/factorial(j)) end do

(1+x)^3

  

3*(1+x)^2*u[1](x)

  

3*(1+x)*u[1](x)^2+3*(1+x)^2*u[2](x)

  

u[1](x)^3+6*(1+x)*u[1](x)*u[2](x)+3*(1+x)^2*u[3](x)

 (12)
> 

A[0]

(1+x)^3

 (13)
> 

y[0] := 1+x

1+x

 (14)
> 

y[1] := invlaplace(2*laplace(A[0], x, s)/s^2, s, x)-invlaplace(x*laplace(diff(y[0], x), x, s)/s^2, s, x)+invlaplace(laplace(y[0], x, s)/s^2, s, x)-invlaplace(laplace(diff(y[0], x), x, s)/s^2, s, x)

(1/10)*x^2*(x^3+5*x^2+10*x+10)-(1/2)*x^3+(1/6)*x^2*(x+3)-(1/2)*x^2

 (15)
> 

y[1] := expand((1/10)*x^2*(x^3+5*x^2+10*x+10)-(1/2)*x^3+(1/6)*x^2*(x+3)-(1/2)*x^2)

(1/10)*x^5+(1/2)*x^4+(2/3)*x^3+x^2

 (16)
> 

"u[1](x) :=y[1] "

proc (x) options operator, arrow, function_assign; y[1] end proc

 (17)
> 

NULL
> 

A[1]

3*(1+x)^2*((1/10)*x^5+(1/2)*x^4+(2/3)*x^3+x^2)

 (18)
> 

y[2] := invlaplace(2*laplace(A[1], x, s)/s^2, s, x)-invlaplace(x*laplace(diff(y[1], x), x, s)/s^2, s, x)+invlaplace(laplace(y[1], x, s)/s^2, s, x)-invlaplace(laplace(diff(y[1], x), x, s)/s^2, s, x)

(1/840)*x^4*(7*x^5+63*x^4+212*x^3+476*x^2+672*x+420)-(1/60)*x^4*(x^3+6*x^2+10*x+20)+(1/420)*x^4*(x^3+7*x^2+14*x+35)-(1/60)*x^3*(x^3+6*x^2+10*x+20)

 (19)
> 

y[2] := expand(%)

(1/120)*x^9+(3/40)*x^8+(5/21)*x^7+(7/15)*x^6+(17/30)*x^5+(1/12)*x^4-(1/3)*x^3

 (20)
> 

" u[2](x):=y[2]"

proc (x) options operator, arrow, function_assign; y[2] end proc

 (21)
> 

A[2]

3*(1+x)*((1/10)*x^5+(1/2)*x^4+(2/3)*x^3+x^2)^2+3*(1+x)^2*((1/120)*x^9+(3/40)*x^8+(5/21)*x^7+(7/15)*x^6+(17/30)*x^5+(1/12)*x^4-(1/3)*x^3)

 (22)
> 

y[3] := invlaplace(2*laplace(A[2], x, s)/s^2, s, x)-invlaplace(x*laplace(diff(y[2], x), x, s)/s^2, s, x)+invlaplace(laplace(y[2], x, s)/s^2, s, x)-invlaplace(laplace(diff(y[2], x), x, s)/s^2, s, x)

(1/10810800)*x^5*(7623*x^8+99099*x^7+518778*x^6+1634490*x^5+3647930*x^4+5167305*x^3+4221360*x^2+900900*x-1081080)-(1/25200)*x^5*(21*x^6+210*x^5+750*x^4+1680*x^3+2380*x^2+420*x-2100)+(1/831600)*x^5*(63*x^6+693*x^5+2750*x^4+6930*x^3+11220*x^2+2310*x-13860)-(1/25200)*x^4*(21*x^6+210*x^5+750*x^4+1680*x^3+2380*x^2+420*x-2100)

 (23)
> 

y[3] := expand(y[3])

(286/945)*x^9+(131/336)*x^8+(1/7)*x^10+(17/70)*x^7-(1/40)*x^6-(1/20)*x^5+(1/12)*x^4+(1091/23100)*x^11+(11/15600)*x^13+(11/1200)*x^12

 (24)
> 

NULL
> 

addingterm := y[0]+y[1]+y[2]+y[3]

1+x+(37/60)*x^5+(2/3)*x^4+(1/3)*x^3+x^2+(2351/7560)*x^9+(781/1680)*x^8+(101/210)*x^7+(53/120)*x^6+(1/7)*x^10+(1091/23100)*x^11+(11/15600)*x^13+(11/1200)*x^12

 (25)


 

Download aproximate_and_exact_solution.mw

a table like that

 

How calculate Adomian polynomial by different algo...

Asked by: salim-barzani 1610
ode pde differential-equations
November 04 2024
0  22

there is four formula for calculate them which i know them by name of author the first one is adomian second one is (BiazarShafiofAdomian) which one member of mableprimes write code for me,but i don't know how use for all kind function maybe in future i upload this program for fix this issue, the third one is by zhao which is i think i easy for calculate just  i need someone one to wite the program and do some test for some example i  upload some picture in case for getting algorithm to writting and have some example for testing  so  lets see who can do this algorithm is very usfule when we solve ODE or PDE by LDM, also last method is by taking integral have a good method, in this question this algorithm is zhao which is usfull one

Error integrating a procedure...

Asked by: Gonzalo Garcia 155
integration sum add
November 04 2024
1  2

Hi!

I am using a proceure to conpute de integral of a function by he Simpson's rule. My function is defined from a function and a procedure, but I am getting the error  "Error, (in w) invalid input: hfun2 expects its 1st argument, t, to be of type numeric, but received (1/10)*i+1/20"

As you can see in the attaxhed file, I have tried several ways to compute the integral but always returns the above error. Please, can yo help me?

Thanks

forum.mw

Exploring Oscillators: A Simple Pendulum Worksheet...

by: Paras31 250 Maple 2024
ode education differential-equations
November 03 2024
1  0

Over the last few days, I’ve been creating worksheets on oscillators to support my class’s understanding of these fundamental physics concepts. I wanted to share one of these worksheets that I found particularly useful for illustrating energy exchange and motion dynamics.

A simple pendulum is a classic physics example that exhibits periodic motion. It consists of a mass m (called a bob) attached to a string or rod of length L, which swings back and forth under the influence of gravity. When the bob is displaced from its equilibrium position and released, it swings back and forth under the influence of gravity.
To derive the equation of motion, we can examine the forces acting on the pendulum bob and use Newton’s second law.

Period of a Pendulum:

• 

Frequency (f) "-" the number of cycles the pendulum completes in one second. Measued in hertz ("Hz)."

f = 1/T

• 

Period ("T) -" the time it takes the pendulum to complete one cycle. Measued in seconds (s).

T = 2*Pi*sqrt(L/g)

This period depends only on the length Land gravitational acceleration "g,"meaning it is independent of the amplitude for small oscillations.

What is the period and the frequency of a single pendulum that is 70 cm long on the earth and on the moon?

> 

L := .7; g__earth := 9.8; g__moon := 1.6

.7

  

9.8

  

1.6

 (1)
> 

T__earth := 2*Pi*sqrt(L/g__earth); f__earth := L/T__earth

1.679251909

  

.4168522878

 (2)
> 

T__moon := 2*Pi*sqrt(L/g__moon); f__moon := L/T__moon

4.155936442

  

.1684337597

 (3)

The above image is taken from https://www.researchgate.net/publication/365297210_Scientific_counterfactuals_as_make-believe

1. 

Forces on the Pendulum Bob:

The main forces acting on the bob are:

• 

The gravitational force"`f__g`=mg, "acting vertically downward.
• 

The tension Tauin the string, acting along the string toward the pivot point.
2. 

Components of the Gravitational Force:

Since the pendulum swings in an arc, it’s helpful to resolve the gravitational force into two components:

• 

Radial Component (along the string): This component, "`f__y`=mgcostheta ," is countered by the tension in the string and does not contribute to the pendulum’s motion.
• 

Tangential Component (perpendicular to the string): This component, f__z = -`mgsinθ`(restoring force), acts along the arc of the pendulum’s swing and is responsible for its motion.
3. 

Applying Newton's Second Law
Since the tangential component of the gravitational force causes the pendulum’s motion, we can apply Newton's second law in the tangential direction:``

f__X = ma__x

Substituting for f__x and the tangential acceleration a__xNULL

m*(diff(s(t), t, t)) = -`mgsinθ`

where diff(s(t), t, t) = a__x and a__x = diff(x, t, t)

Now, we want to write everything in terms of θ

s = `Lθ`

we obtain:

diff(theta(t), t, t) = -g*`sinθ`/L

Small-Angle Approximation

For small angles (typically) theta  , the approximation

`≈`(sin(theta), theta)(where theta is in radians) can be used. This simplifies the equation:

diff(theta(t), t, t) = -g*theta/L

This equation is now in the form of a simple harmonic oscillator

diff(theta(t), t, t) = -omega^2*theta

where omega = sqrt(g/L)is the angular frequency of the pendulum.

> 

restart; with(plots); with(DEtools)
> 

L := 1; m := .2; g := 9.8

1

  

.2

  

9.8

 (4)
> 

T := 2*Pi*sqrt(L/g)

2.007089924

 (5)
> 

omega := sqrt(g/L)

3.130495168

 (6)
> 

ODE__1 := diff(theta(t), t, t)+omega^2*theta = 0; IC := theta(0) = A, (D(theta))(0) = 0

diff(diff(theta(t), t), t)+9.799999997*theta(t) = 0

  

theta(0) = A, (D(theta))(0) = 0

 (7)
> 

sol := dsolve({IC, ODE__1}, theta(t))

theta(t) = A*cos((1/100000)*97999999970^(1/2)*t)

 (8)
> 

plot_1 := subs(A = 0.873e-1, sol); plotsresult := plot([rhs(plot_1)], t = 0 .. 2, color = [red])

 
> 

`θ_t` := rhs(subs(A = 0.873e-1, sol)); v_t := diff(`θ_t`, t)

0.873e-1*cos((1/100000)*97999999970^(1/2)*t)

  

-0.8730000000e-6*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

 (9)
> 

T := (1/2)*m*L^2*v_t^2; V := m*g*L*(1-cos(`θ_t`)); H := T+V

0.7468864200e-2*sin((1/100000)*97999999970^(1/2)*t)^2

  

1.96-1.96*cos(0.873e-1*cos((1/100000)*97999999970^(1/2)*t))

  

0.7468864200e-2*sin((1/100000)*97999999970^(1/2)*t)^2+1.96-1.96*cos(0.873e-1*cos((1/100000)*97999999970^(1/2)*t))

 (10)
> 

energy_plot := plot([eval(T), eval(V), eval(H)], t = 0 .. 5, color = [red, blue, green], legend = ["Kinetic Energy", "Potential Energy", "Total Energy"], title = "Energy Exchange in Simple Pendulum", labels = ["Time (s)", "Energy (Joules)"])

 
> 

directionfield := DEplot([diff(theta(t), t) = v(t), diff(v(t), t) = -omega^2*theta(t)], [theta(t), v(t)], t = 0 .. 20, theta = -20 .. 20, v = -40 .. 40, arrows = medium, title = "Direction Field for Simple Harmonic Oscillator", axes = boxed, color = navy)
> 

sol1 := dsolve({ODE__1, theta(0) = 3, (D(theta))(0) = 0}, theta(t)); sol2 := dsolve({ODE__1, theta(0) = 6.5, (D(theta))(0) = 0}, theta(t)); sol3 := dsolve({ODE__1, theta(0) = -8, (D(theta))(0) = 0}, theta(t)); sol4 := dsolve({ODE__1, theta(0) = 9.7, (D(theta))(0) = 2.5}, theta(t))

theta(t) = 3*cos((1/100000)*97999999970^(1/2)*t)

  

theta(t) = (13/2)*cos((1/100000)*97999999970^(1/2)*t)

  

theta(t) = -8*cos((1/100000)*97999999970^(1/2)*t)

  

theta(t) = (25000/9799999997)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)+(97/10)*cos((1/100000)*97999999970^(1/2)*t)

 (11)
> 

theta1 := rhs(sol1); theta2 := rhs(sol2); theta3 := rhs(sol3); theta4 := rhs(sol4)

3*cos((1/100000)*97999999970^(1/2)*t)

  

(13/2)*cos((1/100000)*97999999970^(1/2)*t)

  

-8*cos((1/100000)*97999999970^(1/2)*t)

  

(25000/9799999997)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)+(97/10)*cos((1/100000)*97999999970^(1/2)*t)

 (12)
> 

v1 := diff(theta1, t); v2 := diff(theta2, t); v3 := diff(theta3, t); v4 := diff(theta4, t)

-(3/100000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

  

-(13/200000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

  

(1/12500)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

  

(5/2)*cos((1/100000)*97999999970^(1/2)*t)-(97/1000000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

 (13)
> 

phase_plot := plot([[eval(theta1, t = tval), eval(v1, t = tval), tval = 0 .. 20], [eval(theta2, t = tval), eval(v2, t = tval), tval = 0 .. 20], [eval(theta3, t = tval), eval(v3, t = tval), tval = 0 .. 20], [eval(theta4, t = tval), eval(v4, t = tval), tval = 0 .. 20]], style = line, title = "Phase Portrait for Simple Harmonic Oscillator", labels = ["x (Displacement)", "v (Velocity)"], color = ["red", "blue", "green", "orange"], axes = boxed)
> 

display(directionfield, phase_plot)

 
> 

theta_at_1_sec := subs(t = 1, A = 0.873e-1, rhs(sol)); evalf(theta_at_1_sec)

0.873e-1*cos((1/100000)*97999999970^(1/2))

  

-0.8729462437e-1

 (14)
 

NULL

Download Simple_Pendulum.mw

How to solve this equation?...

Asked by: mary120 75
November 03 2024
0  3

Hi,
I have an equation and I want to solve it parametrically to find x , but I couldn't do that with "solve" command. (I know x should be  real and positive). What should I do?
Root_of.mw

How collect series solution of PDE?...

Asked by: salim-barzani 1610
pde differential-equations laplace
November 02 2024
1  31

this example is easiest one for getting solution but i can't collect each part and do like elite i can do each part seperatly but it take to much time i want collect solution and do by easier way if possible this is a laplace adomian decomposition methd which contain adomian polynomial too i want upgrade the code, can any one help the  process for get better vision of this topic 
i do upload some picture and my mw. for more undrestanding

and please can any one explan why when i take laplace why is write D[2](u)(x, 0) must be D[1]?

> 

restart
> 

with(inttrans)
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(SolveTools)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

 (1)
> 

with(PDEtools)
> 

declare()

`Declared :`

  

u(x, t)*`to be displayed as`*u

  

`The prime differentiation variable has not been declared yet`

  

`Displayed derivatives and declared functions will be copied and pasted "as they were entered"`

 (2)
> 

declare(u(x, t))

u(x, t)*`will now be displayed as`*u

 (3)
> 

eq := diff(u(x, t), `$`(t, 2))+u(x, t)^2-(diff(u(x, t), x))^2 = 0

diff(diff(u(x, t), t), t)+u(x, t)^2-(diff(u(x, t), x))^2 = 0

 (4)
> 

NULL
> 

eqs := laplace(eq, t, s)

s^2*laplace(u(x, t), t, s)-(D[2](u))(x, 0)-s*u(x, 0)+laplace(u(x, t)^2, t, s)-laplace((diff(u(x, t), x))^2, t, s) = 0

 (5)
> 

solve({eqs}, {laplace(u(x, t), t, s)})

{laplace(u(x, t), t, s) = (s*u(x, 0)+(D[2](u))(x, 0)-laplace(u(x, t)^2, t, s)+laplace((diff(u(x, t), x))^2, t, s))/s^2}

 (6)
> 

subs({u(x, 0) = 0, (D[2](u))(x, 0) = exp(x)}, %)

{laplace(u(x, t), t, s) = (exp(x)-laplace(u(x, t)^2, t, s)+laplace((diff(u(x, t), x))^2, t, s))/s^2}

 (7)
> 

eq3 := invlaplace(%, s, t)

{u(x, t) = exp(x)*t-(int(u(x, _U1)^2*(t-_U1), _U1 = 0 .. t))+int((diff(u(x, _U1), x))^2*(t-_U1), _U1 = 0 .. t)}

 (8)
> 

NULL
> 

NULL
> 

NULL
> 

"u[0](x,t):=exp(x)*t"

proc (x, t) options operator, arrow, function_assign; exp(x)*t end proc

 (9)
> 

n := N

N

 (10)
> 

k := K

K

 (11)
> 

f := proc (u) options operator, arrow; u^2 end proc

proc (u) options operator, arrow; u^2 end proc

 (12)
> 

for j from 0 to 3 do A[j] := subs(lambda = 0, (diff(f(seq(sum(lambda^i*u[i](x, t), i = 0 .. 20), m = 1 .. 2)), [`$`(lambda, j)]))/factorial(j)) end do

(exp(x))^2*t^2

  

2*exp(x)*t*u[1](x, t)

  

u[1](x, t)^2+2*exp(x)*t*u[2](x, t)

  

2*u[1](x, t)*u[2](x, t)+2*exp(x)*t*u[3](x, t)

 (13)
> 

NULL
> 

NULL
> 

n := N

N

 (14)
> 

k := K

K

 (15)
> 

f := proc (u) options operator, arrow; (diff(u, x))^2 end proc

proc (u) options operator, arrow; (diff(u, x))^2 end proc

 (16)
> 

for j from 0 to 3 do B[j] := subs(lambda = 0, (diff(f(seq(sum(lambda^i*u[i](x, t), i = 0 .. 20), m = 1 .. 2)), [`$`(lambda, j)]))/factorial(j)) end do

(exp(x))^2*t^2

  

2*exp(x)*t*(diff(u[1](x, t), x))

  

(diff(u[1](x, t), x))^2+2*exp(x)*t*(diff(u[2](x, t), x))

  

2*(diff(u[1](x, t), x))*(diff(u[2](x, t), x))+2*exp(x)*t*(diff(u[3](x, t), x))

 (17)
> 

"#` know we need find all term of u[0]=exp(x)*t` #` u`[1]=-invlaplace(1/(s^(2))(`A__0`+`B__0`)) u[2]=-invlaplace(1/(s^(2))(`A__1`+`B__1`)) ans so on ... at end i want collect all of them and find final result even if is aproximate and want do test of pde too "
> 

NULL

Download explananing_of_get_solution.mw

Typesetting internal error when Physics update is ...

Asked by: nm 11398
libname physics latex
November 02 2024
0  3

[moderator: The Physics update Library fixes this bug with the same error generated and reported by the same Mapleprimes member on another ODESteps example.]

I have removed Physics update from libname path. 

Now I find I get error calling latex command. When Physics update is on libname, no error.

The error is 

         Error, (in Typesetting:-Parse) too many levels of recursion

I am using worksheet with typesetting extended. But also when I change it to typesetting standard, same error. 

Does this mean one must keep Physics update on libname path for Maple to work OK?

Is this error expected if Physics update is not on libname?

Worksheet below that shows this problem

> 

restart;
> 

interface(version);

`Standard Worksheet Interface, Maple 2024.2, Windows 10, October 29 2024 Build ID 1872373`

CASE 1. With Physics update in libname path, no error

 
> 

restart;

> 

libname;

"C:\Users\Owner\maple\toolbox\2024\Physics Updates\lib", "C:\Program Files\Maple 2024\lib"

> 

ode:=[diff(x__1(t),t)=2*x__1(t)+x__2(t),diff(x__2(t),t)=2*x__1(t)+3*x__2(t)];

[diff(x__1(t), t) = 2*x__1(t)+x__2(t), diff(x__2(t), t) = 2*x__1(t)+3*x__2(t)]

> 

the_output:=Student:-ODEs:-ODESteps(ode,output=typeset):

> 

latex(the_output,'output'=string):

> 

 

CASE 2.  Removing Physics from libname path, gives internal error

 
> 

restart;

> 

libname:=libname[2];

"C:\Program Files\Maple 2024\lib"

> 

ode:=[diff(x__1(t),t)=2*x__1(t)+x__2(t),diff(x__2(t),t)=2*x__1(t)+3*x__2(t)];

[diff(x__1(t), t) = 2*x__1(t)+x__2(t), diff(x__2(t), t) = 2*x__1(t)+3*x__2(t)]

> 

the_output:=Student:-ODEs:-ODESteps(ode,output=typeset):

> 

latex(the_output,'output'=string):

Error, (in Typesetting:-Parse) too many levels of recursion

> 

 

 

 

Download internal_error_from_latex_when_libname_changed_nov_2_2024.mw

plotting numerically solved PDE at specific range ...

Asked by: alfarunner 90
pde numeric pdsolve
November 02 2024
0  1

it's easy to plot a numeric solution of a pde, e.g u(x,t) at a specific time but it's more difficult to plot a solution for a specific x-value as a function of time. How do I do that in Maple?

I find an error...

Asked by: Yuri Zamorano 30
regression statistics incomplete-question
November 01 2024
0  1

I find an error when calculating ExponentialFit(X, Y, x, summarize = embed) function. It's have an error in R-squared and Adjusted R-squared.

Issues with Variational Derivative and Commutator...

Asked by: WA573 80
November 01 2024
0  4

I’m working on deriving the equation of motion (16) from equation (9) (see attached image), but I’ve encountered a couple of issues in Maple. i) The variational derivative of Lagrangian, and ii) the commutator is zero, which isn’t expected.

How do we fix these issues? Are there any specific packages regarding the construction of Lagrangian?

VariationDerivative.mw

How to compute with "exports" from dsolve,numeric...

Asked by: C_R 3442
ode procedure dsolve
November 01 2024
1  12

As an example, I have tried to compute the length of a pendulum to verify how good the lagrange multiplier establishes the constraint of a constant length.

Computation of single values from a formula was possible. Using the same formulas as argument of Maple commands (plot in this case) turned out to be complicated.

I tried many variants with unevaluation quotes before turning to a call to unapply within unevaluation quotes.
There must be simpler ways. Please advise.

PS.: I have put exports in quotes because I assume it would have been easier if dsolve,numeric could export functions. Am I correct in this assumption?

From dsolve,numeric,DAE

> 

dsys := {1 - x(t)^2 - y(t)^2, diff(x(t), t, t) = -2*x(t)*lambda(t), diff(y(t), t, t) = -9.8 - 2*y(t)*lambda(t)}

{1-x(t)^2-y(t)^2, diff(diff(x(t), t), t) = -2*x(t)*lambda(t), diff(diff(y(t), t), t) = -9.8-2*y(t)*lambda(t)}

 (1)
> 

dini := {lambda(0) = 5.025000000, x(0) = 0, y(0) = -1, D(x)(0) = 1/2, D(y)(0) = 0}

{lambda(0) = 5.025000000, x(0) = 0, y(0) = -1, (D(x))(0) = 1/2, (D(y))(0) = 0}

 (2)
> 

dsol1 := dsolve(dsys union dini, numeric);

proc (x_rkf45_dae) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45_dae) else _xout := evalf(x_rkf45_dae) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 28, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..65, {(1) = 5, (2) = 4, (3) = 0, (4) = 0, (5) = 0, (6) = 2, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.5047658755841546e-2, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..5, {(1) = .0, (2) = .5, (3) = -1.0, (4) = .0, (5) = 5.025}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..5, {(1) = .1, (2) = .1, (3) = .1, (4) = .1, (5) = .1}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = -2.0, (1, 4) = .0, (2, 1) = .5, (2, 2) = .0, (2, 3) = .0, (2, 4) = -1.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..4, 1..4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0}, datatype = float[8], order = C_order), Array(1..5, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8]), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..10, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = 0, (2) = 0, (3) = 0, (4) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..5, {(1) = .0, (2) = .5, (3) = -1.0, (4) = .0, (5) = undefined}, datatype = float[8], order = C_order), Array(1..5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .5, (2) = -.0, (3) = .0, (4) = .25}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) local Z1, Z2, Z3, Z4, Z5; Z1 := Y[2]^2; Z2 := Y[4]^2; Z3 := Z1+Z2; Z4 := Y[1]^2; Z5 := Y[3]^2+Z4; Z5 := 1/Z5; Y[5] := -(1/10)*(-5*Z3+49*Y[3])*Z5; if N < 1 then return 0 end if; YP[2] := -(1/5)*(5*Z3-49*Y[3])*Y[1]*Z5; YP[4] := -(1/5)*(5*Y[3]*(Z1+Z2)+49*Z4)*Z5; YP[1] := Y[2]; YP[3] := Y[4]; 0 end proc, -1, proc (X, Y, R) R[1] := Y[1]^2+Y[3]^2-1; R[2] := Y[1]*Y[2]+Y[3]*Y[4]; 0 end proc, proc (X, Y, J) J[1 .. 2, 1 .. 4] := 0; J[1, 1] := 2*Y[1]; J[1, 3] := 2*Y[3]; J[2, 1] := Y[2]; J[2, 2] := Y[1]; J[2, 3] := Y[4]; J[2, 4] := Y[3]; 0 end proc, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([proc (X, Y, R) R[1] := Y[1]^2+Y[3]^2-1; R[2] := Y[1]*Y[2]+Y[3]*Y[4]; 0 end proc, proc (X, Y, J) J[1 .. 2, 1 .. 4] := 0; J[1, 1] := 2*Y[1]; J[1, 3] := 2*Y[3]; J[2, 1] := Y[2]; J[2, 2] := Y[1]; J[2, 3] := Y[4]; J[2, 4] := Y[3]; 0 end proc, 0, 0, Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), [-1+x(t)^2+y(t)^2, x(t)*(diff(x(t), t))+y(t)*(diff(y(t), t))]]), ( 17 ) = ([proc (N, X, Y, YP) local Z1, Z2, Z3, Z4, Z5; Z1 := Y[2]^2; Z2 := Y[4]^2; Z3 := Z1+Z2; Z4 := Y[1]^2; Z5 := Y[3]^2+Z4; Z5 := 1/Z5; Y[5] := -(1/10)*(-5*Z3+49*Y[3])*Z5; if N < 1 then return 0 end if; YP[2] := -(1/5)*(5*Z3-49*Y[3])*Y[1]*Z5; YP[4] := -(1/5)*(5*Y[3]*(Z1+Z2)+49*Z4)*Z5; YP[1] := Y[2]; YP[3] := Y[4]; 0 end proc, -1, proc (X, Y, R) R[1] := Y[1]^2+Y[3]^2-1; R[2] := Y[1]*Y[2]+Y[3]*Y[4]; 0 end proc, proc (X, Y, J) J[1 .. 2, 1 .. 4] := 0; J[1, 1] := 2*Y[1]; J[1, 3] := 2*Y[3]; J[2, 1] := Y[2]; J[2, 2] := Y[1]; J[2, 3] := Y[4]; J[2, 4] := Y[3]; 0 end proc, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0), ( 28 ) = (0) ] )) ] ); _y0 := Array(0..5, {(1) = 0., (2) = 0., (3) = .5000000000000000000000, (4) = -1., (5) = 0.}); _vmap := array( 1 .. 5, [( 1 ) = (5), ( 2 ) = (1), ( 3 ) = (2), ( 4 ) = (3), ( 5 ) = (4) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, lambda(t), x(t), diff(x(t), t), y(t), diff(y(t), t)], (4) = []}); _vars := _dat[3]; _pars := map(lhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45_dae, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45_dae, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45_dae, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45_dae, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45_dae), 'string') = rhs(x_rkf45_dae); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45_dae), 'string') = rhs(x_rkf45_dae)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45_dae) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45_dae) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (3)
> 

dsol1(1/2);

[t = .500000000000000, lambda(t) = HFloat(4.837512079685444), x(t) = HFloat(0.15920394558985643), diff(x(t), t) = HFloat(0.003966955001051789), y(t) = HFloat(-0.9872457162788973), diff(y(t), t) = HFloat(6.396972656506275e-4)]

 (4)
> 

r:=sqrt(x(t)^2+y(t)^2);#expression to compute

(x(t)^2+y(t)^2)^(1/2)

 (5)
> 

plots:-odeplot(dsol1,[t,r],t=0..10)

 
> 

eval(r,dsol1(0.5));

HFloat(1.0000000003012055)

 (6)
> 

plot(eval(r,dsol1(t)),t=0..10)

Error, invalid input: eval received dsol1(t), which is not valid for its 2nd argument, eqns

  
> 

plot('eval(r,dsol1(t))',0..10)

Error, (in plot) procedure expected, as range contains no plotting variable

  
> 

plot(unapply('eval(r,dsol1(t))',t),0..10);#complicated way

 
> 

r:=t->sqrt(x(t)^2+y(t)^2);#function to compute

proc (t) options operator, arrow; sqrt(x(t)^2+y(t)^2) end proc

 (7)
> 

f:=x-> dsol1(x),'r(t)'

proc (x) options operator, arrow; dsol1(x) end proc, r(t)

 (8)
> 

subs[eval](f(0.5))

HFloat(1.0000000003012055)

 (9)
> 

plot(unapply('subs[eval](f(t))',t),0..10);# equally complicated

 

NULL

Download computation_with_dsolve_exports.mw

Edit: I found these after sending but I still think that there is more to improve. Is eval aor subs needed at all?

plot('eval(r,dsol1(t))',t=0..10)
plot(''eval(r,dsol1(t))'',t=0..10)
First 62 63 64 65 66 67 68 Last Page 64 of 2219