Preben Alsholm

MaplePrimes Activity

These are answers submitted by Preben Alsholm

Isn't this just a case of garbage in, garbage out?
Using allvalues gives me 2 matrices:

x := RootOf(z^2-t, z);
m := Matrix(2, (i,j) -> add(a[j, k]*((-1)^(i-1))^k*x^k, k = 0..1));
L1,V1 := Eigenvectors(m1);
## and
L2,V2 := Eigenvectors(m2);


The two outer ifs have no else, therefore if not H[i] < 2.7 or if not A[i, f] < 12 then nothing is going to happen.

So you must tell us (or yourself) what you want to do in those cases.
For clarity I wrote the code with some indentation:

for i to n do 
   if H[i] < 2.7 then 
     if A[i, f] < 12 then 
       if A[i, o] < 1.2 then 
         Q[i, foo] := evalf(610*(A[i, o]*sqrt(H[i])*h[k]*A[i, T])^(1/2)) 
         Q[i, foo] := evalf(7.8*A[i, t]+378*A[i, o]*sqrt(H[i])) 
       end if 
     end if 
   end if; 
   print(Q[i, foo]) 
end do;


You seem to have some syntactical problems with your derivatives. Being no friend of 2D math input I made use of Kitonum's translation into the trustworthy 1D math input, known also as Maple Input.

interface(imaginaryunit = j);
k := .5; tau := .95; mu := 0.1e-1; pi := 116.1; theta := 0.8e-2; phi := 0.25e-2; epsilon := 0.2e-2; rho := 0.5e-1; beta := 0.115e-1; chi := 0.598e-2; q := .5; eta := .2; delta := .1; alpha := 0.57e-1; p := .2; Upsilon := 1.2;
lambda := k*tau*(C(t)*Upsilon+(I)(t))/(S(t)+V(t)+C(t)+(I)(t)+R(t));
sys:={diff(C(t), t) = rho*lambda*S(t)+rho*epsilon*lambda*V(t)+(1-q)*eta*I(t)-(mu+beta+chi)*C(t), diff(I(t), t) = (1-rho)*lambda*S(t)+(1-rho)*epsilon*lambda*V(t)+chi*C(t)-(mu+alpha+eta)*I(t), diff(R(t), t) = beta*C(t)+q*eta*I(t)-(mu+delta)*R(t), diff(S(t), t) = (1-p)*pi+phi*V(t)+delta*R(t)-(mu+lambda+theta)*S(t), diff(V(t), t) = p*pi+theta*S(t)-(epsilon*lambda+mu+phi)*V(t)};
ics:= {S(0) = 8200, V(0) = 2800, C(0) = 200, (I)(0) = 210, R(0) = 200};
res:=dsolve(sys union ics, numeric); 

The 3 dimensional plot of C,I,S is shown here. Notice that I only plot after the initial settling has taken place. I chose 2000..3000:

Below we see that your system has 3 equilibrium points and that they all are unstable.

for pp in pts do simplify(LinearAlgebra:-Eigenvalues(J(op(pp)))) end do;


In maple z^(1/3) is the principal cube root of z. Thus as an example (-1)^(1/3) is equal to 1/2+I*sqrt(3)/2.
Try evalc( (-1)^(1/3) );
What you can use since you expect to get -1 from (-1)^(1/3) is surd.

Thus your plot will be produced by


Your code has some syntactical problems. I tried to clear them up.
Besides that the second derivatives w.r.t. eta cannot appear in the boundary conditions. I removed them.
Thus the code looks like this:

Digits := 10;
sys_ode := 2*diff(T(eta), eta, eta, eta)+T(eta)*diff(T(eta), eta, eta) = 0;
ics := T(0) = 0, D(T)(0) = 1, D(T)(20) = 0:
sol1 := dsolve({ics, sys_ode}, numeric, output = operator);
PDE1 := eval([2/P * diff(g(x, eta), eta, eta)+q(eta)*diff(g(x, eta), eta)-g(x, eta)*w(eta) = 2*x*w(eta)*diff(g(x, eta), x), 2/S * diff(phi(x, eta), eta, eta)+q(eta)*diff(phi(x, eta), eta) = 2*x*w(eta)*diff(phi(x, eta), x)]);
subBC1 := -phi(x, 0)*exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)));
subBC2 := alpha*phi(x, 0)*sqrt(x)*exp(g(x,0)*sqrt(x)/(1+varepsilon*g(x,0)*sqrt(x)));
BC := {g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, D[2](g)(x, 0) = subBC1, D[2](phi)(x, 0) = subBC2};
P := 1;
S := 1;
varepsilon:=1: #Added
alpha:=1: #Added
pds := pdsolve(PDE1, BC, numeric, spacestep = .25);

I receive the error:

Error, (in pdsolve/numeric/match_PDEs_BCs) unable to associate boundary conditions and dependent variables, system is too complex
Even when setting varepsilon=0 the same error comes up.
In an test to see if there may be other problems I tried these boundary conditions, where exp is replaced by 1:

BCtest:={g(0, eta) = 0, g(x, 20) = 0, phi(0, eta) = 1, phi(x, 20) = 1, D[2](g)(x, 0) = -phi(x, 0), D[2](phi)(x, 0) = alpha*sqrt(x)*phi(x, 0)};

The answer came up right away with no complaint, so it appears that the complexity of BC is the problem.

You could try simplify/size:

length(Teller); # 4308
length(TS); # 2035
whattype(TS); # `+`
nops(TS); # 5

Begin by not assigning to the parameters you want to change later.

lambda := 1.0; K__r := 1.0; S__r := .5; m := .5; M := sqrt(10.0); `&varkappa;` := .5; Omega := sqrt(5.0); 
## Have omitted Gm, Gr, P__r, S__c
## Now define PDE and IBC as you did
##  ...
## Define SOL as a procedure:
SOL:=(gm,gr,pr,sc)->pdsolve(eval(PDE,{Gm=gm,Gr=gr,P__r=pr,S__c=sc}), IBC, numeric, spacestep = 0.1e-1);
## Test:
SOL(5,6,0.71,0.22):-plot(t = .3, color = red);
SOL(5,6,0.71,0.1):-plot(t = 1,color = red);

Adding this code you can get an animation;

PL1:=sc->SOL(5,6,0.71,sc):-plot(t = 1,color = red);

I'm assuming that you meant sys:= {...} and not sys*{...}.
Thus we have:

sys:={A2+D2 = 0, B1*sin(192*K1)+D1 = 0, 3.383*B2*K1+C2 = 0, B1*K1^2*sin(192*K1)+11.444689*A2*K1^2*cos(568.344*K1)+11.444689*B2*K1^2*sin(568.344*K1) = 0, A2*cos(568.344*K1)+B2*sin(568.344*K1)+168*C2+D2 = 0, -3.383*A2*K1*sin(568.344*K1)+3.383*B2*K1*cos(568.344*K1)+C2-B1*K1*cos(192*K1) = 0};
##You cannot use fsolve if you want to keep K1 unassigned.
solve(sys,{A2, B1, B2, C2, D1, D2});

You just get the trivial solution {A2 = 0., B1 = 0., B2 = 0., C2 = 0., D1 = 0., D2 = 0.}.
## You may ask: For which values of K1 are there nontrivial solutions?
Notice that your equations are linear in A2, B1, B2, C2, D1, D2.
To look into that use a little linear algebra:

A,zz:=LinearAlgebra:-GenerateMatrix(sys,[A2, B1, B2, C2, D1, D2]);
dt:=LinearAlgebra:-Determinant(A); #Must be zero to get nontrivial solutions
fsolve(dt=0,K1=0); # Answer 0.
fsolve(dt=0,K1=0.006); # Answer: 0.5941860004e-2

We must expect infinitely many values of K1 with that property.
# Since dt-eval(dt,K1=-K1) evaluates to zero, you need only look at K1 >= 0.
To find the corresponding solutions for A2, B1, B2, C2, D1, D2 you could do:

K1a:=fsolve(dt=0,K1=0.006); #Example only
G0:=fnormal~(G); #When dt=0 the rank < full. Removing roundoff.



I changed pi to Pi.

This probably won't be satisfactory, but I tried without using plotsetup. I right clicked on the plot and exported as a jpg file:

plots:-logplot(x^2-3, x = 0 .. 100, font = [bold, "TimesNewRoman", 30],size=[1500,1100]);

The plot:

The functions named _F1, _F2, ... are meant as arbitrary functions. Thus an expression like _F1(t+x) means any expression depending on the sum t+x only. Nothing more is meant (except that appropriate smoothness requirements are implicitly assumed).

You need initial conditions for the derivative too. Your ode is called ode1, not as you have ode.
To get arrows you must rewrite as a system of first order:

ode1 := diff(x(t), t, t) = (5*9.80665)*sin((1/6)*Pi)-(10*(10-sqrt(x(t)^2+25)))*x(t)/sqrt(x(t)^2+25)-(diff(x(t), t));
ics:=[seq([x(0) = (1/4)*k,D(x)(0)=k/2], k = -20 .. 20)]; # I chose D(x)(0)=k/2
DEtools[DEplot](ode1, x(t), t = -2 .. 2, ics, x=-7..20, stepsize = 0.5e-1, linecolour = red);
sys := {ode2,subs(ode2,ode1)};
DEtools[DEplot](sys, [x(t),y(t)], t = -2 .. 2, ics_xy, x=-7..15,y=-7..20,color = blue, stepsize = 0.5e-1, linecolour = red,arrows=medium);

The last plot:

As we all agree this is just a "cleaning up" business. There is nothing wrong with the solution.
I just tried to make a procedure CleanUp. It assumes that the input is a solution to a linear ode and comes from dsolve.
It hasn't been tested thoroughly.
For what it is worth here it is:

CleanUp:=proc(sol::equation) local t,cs,Hom,Inhom,S,tm;
  if not lhs(sol)::function(name) then error "This is not an output from dsolve" end if;
  if cs={} then return sol end if;
  if not type(Inhom,`+`) then return sol end if;
  for tm in Inhom do
    if solve('identity'(Hom=tm,t))<>NULL then S:=S union {tm} end if
  end do;
end proc:


You seem to have answered your own question.
Try this:

Student:-NumericalAnalysis:-Quadrature(sin(x), x=0..Pi, method=newtoncotes[open,4], output=information);

There is a discrepancy, however, between what is written in the help page
under Options and what is actually the case.
We read:
"adaptive = true or false
  Whether to apply the adaptive version of the specified method. The following methods are available when adaptive is set to true:
Newton-Cotes Rule, open or closed, with a user-specified order"

but if you try

Student:-NumericalAnalysis:-Quadrature(sin(x), x=0..Pi, adaptive=true,method=newtoncotes[open,4], output=information);

then you will get the error message

Error, (in Student:-NumericalAnalysis:-Quadrature) the newtoncotes[open, 4] method is not available when the adaptive option is specified

Note: I have submitted an SCR.



I changed pi into Pi and made an assignment to K in your worksheet.
I chose as an example to make all parameters in the system equal to 1 and all initial values equal to zero.
Change that as appropriate.
To get the plots of all 4 separately you can use this code:


The plot at the end showing V1 and V2 only (and together):

First 12 13 14 15 16 17 18 Last Page 14 of 149