## 1400 Reputation

17 years, 294 days
University of Twente (retired)
Enschede, Netherlands

My "website" consists of a Maple Manual in Dutch

## s := s+1?...

You seem to use a counter s in the inner loop? Why? It is not initialized, so it causes a recursive assignment error message.
I commented out some unnecessary statements:

`for k to 22 do   T[k] := .8*T[k-1];   for i to n do f(x[i-1]);     # f[i-1];     y[i] := x[i-1]+b[i];     # f(y[i]);     if f(y[i]) < f(x[i-1]) then       x[i] := y[i];  # s := s+1     elif u[i] <= exp(-(f(y[i])-f(x[i-1]))/T[k-1]) then x[i] := y[i]     else x[i] := x[i-1]     end if;     x[i]   end do end do;`

## The function a(u)...

You mean a local maximum of a as a function of u? You can easily solve for a:

`restart;l:=1:alpha:=1:b:=100:k:=20:eq := (alpha+(l+alpha)*u+alpha*k*u^2)*a =       u*(alpha+(l+alpha)*u+alpha*k*u^2)*(1+l*alpha*b/((alpha+(l+alpha)*u+alpha*k*u^2))):A := unapply( solve( eq, a ), u );plot( A, 0..12 ); # There is a local maximum for u < 1U := fsolve( D(A)(u), u=0..1 ):U,A(U);                   0.2306248651, 9.364260108`

## dsolve...

See ?dsolve

`restart;DE := diff(v(t),t) - v(t)/(15-t) = 300/(15-t) - g;s := dsolve( {DE,v(0)=0}, v(t) );subs( g=9.81, s ): V := unapply( subs( %, v(t) ), t );plot( V, 0..14, view=0..3200 );`

## Initialization...

Main problem is that you don't initialize x1 and x2. This causes a "recursive assignment" error.

`restart:N := 10: h := 0.1: X := Array(1..N,0..1): x1 := 1: x2 := 3:for i from 1 to N do t:=h*i: X[i,0]:=x1+h*(-x1)*x2+x1: X[i,1]:=x2+h*x1*x2-x2: x1:=X[i,0]: x2:=X[i,1]: print('Iteration',i,'AtTime',t,x1,x2)end do:X;`

Do you really want X to be an Array? A Matrix (indices 1..N, 1..2) gives a far more readable output:

`convert(X,Matrix);`

(1) You use t and tau for the independent variable; choose one of these

(2) The same for the dependent variable q and Q. Are these the same functions?

(3) Syntax errors, for instance

diff(q(t), t) = ax*(-y(t)-z(t)- ...

Perhaps you mean:

diff(q(t), t) = a*x(t)*(-y(t)-z(t)-

(4) There are several unassigned parameters (a, c, beta, eps_y, ...). These must be assigned for applying dsolve,numeric

(5) And the most serious problem: I'm afraid that nested functions like x(Q(tau)) don't work...

## Remember table...

You can make a function by only defining a remember table:

`for i to 64 do f(A[i]) := B[i] end do:eval(f);     proc () option remember,  'procname(args)' end proc`

## Your original question...

If you differentiate the unevaluated integral Int(q(i),i=0..M) w.r.t. q(i), you get Int(1, i=0..M). This can be achieved by using q instead of q(i).
But later you asked:

"Originally, I wanted to calculate this.

v:=Int(q(i)^2,i=0..M)^2;
dv:=diff(v,q(i));

So you can do:

`v:=Int(q^2,i=0..M)^2;diff(v,q); subs( q=q(i), % ): dv := simplify(%);`

## What do you mean?...

perhaps

`M1 := <<12, 12.5, 13, 13.5, 14, 14.5, 15>|<1.622712644, 1.265443137, 1.028604736,                 .8605013333, .7352916667, .6386248233, .5618945274>>:M2 := <<12, 12.5, 13, 13.5, 14, 14.5, 15>|<5.483608580, 4.289400489, 3.496793877,                2.933480578, 2.513320599, 2.188469637, 1.930230220>>:plot([M1,M2]);`

If not, make clear what you really want.

## Some hints...

(1) You use s and S. Note that small and capital S are different variables

(2) Don't use capital I as a variable name, because that is the imaginary unit

(3) S, L, I and T are the functions that you want to solve for, so the first equation may become:

diff(S(t),t) = Lambda_0 - beta*c*S(t)*J(t)/N - mu*S(t)

(4) You need initial conditions.

## Lookup table...

It is not quite clear what you mean. What if y is in another order? If you want the y's in the order indicated by x, you may construct a lookup table:

`x := [2,3,5,1,4]:y := [1,2,3,4,5]:for i to nops(x) do H[x[i]] := y[i] end do:convert(H,list);                        [4, 1, 2, 5, 3]restart;x := [2,3,5,1,4]:y := [1,2,5,4,3]:for i to nops(x) do H[x[i]] := y[i] end do:convert(H,list);                        [4, 1, 2, 3, 5]`

## Another model...

Carl Love's model doesn't give a good fit, to put it mildly:

`f := a+b*x^c:q := Statistics:-NonlinearFit(f, X, Y, x):display( {pointplot(Data),plot(q,x=50..2000)} );`

Seems to be a serious bug in NonlinearFit.
With the LestSquares in the CurveFitting package, we can make the function a*exp(b*x), by taking logarithms:

` LL := x -> [x[1],ln(x[2])]:nData := LL~(Data):q := CurveFitting:-LeastSquares( nData, x ):Q := exp(op(1,q))*exp(op(2,q));display( {pointplot(Data),plot(Q,x=50..2000)} );`

## Try this...

Do you mean something like

`restart;eta := (x,y) -> y*sqrt(U*v/x);Phi := (x,y) -> sqrt(U*v*x)*f(eta(x,y));e3 := vx = diff(Phi(x,y),x);e4 := vy = diff(Phi(x,y),y);`

## Output=listprocedure...

"But I still would like to know how to differentiate the function r(t) which is numeric"

With the option output=listprocedure, you will get the solution in a format that can be used in a subs-command. Note that it contains not only a procedure for r(t), but also for D(r)(t):

`restart;epsilon := 0.167e-1: AE := 149597870700: mu := 1: v_T := 30290: r_0 := AE: GM := 132712440018*10^9: L_0 := sqrt(GM*mu^2*r_0): d := 86400: t0 := 365.25: M0 := 10^10: M_0 := M0*Heaviside(t0-t): M_1 := M0*Heaviside(t-t0):ode := diff(r(t), t, t)+(-(M_0*t+M_1*t0+L_0)^2/(mu^2*r(t)^3)+GM/r(t)^2)*d^2 = 0:sol := dsolve({ode, r(0) = r_0, (D(r))(0) = 0}, r(t), numeric, output=listprocedure);R := subs( sol, r(t) ): RP := subs( sol, diff(r(t),t) ):`

Now you have the functions R, and its derivative RP:

```plot( RP, 0..500 );
plot(R/AE, 0 .. 500);
```

"Besides: how do I eval r(t) like eval(r(t),t=0) for example."

Simply by R(0);

## Lists, sets and vectors...

There is an important difference between lists (square brackets), sets (curly brackets and sequences (parentheses). See ?list

The input of igcd must be the sequence of the elements of a list, that can be obtained by op(..).
So, for example:

`lst := [[12,8],[16,3]]:(igcd@op)~(lst);                             [4, 1]`

The ?tilde operator works nearly the same as map

## Maxfun and degrees...

(1) See ?dsolve,numeric,maxfun. I can not make it more clear.

(2) You may consider the possibility to make your plot in radians, but to put ?tickmarks in multiples of say 15 degrees along the axis.

 First 7 8 9 10 11 12 13 Last Page 9 of 27
﻿