Advanced Search

Dkunb

65 Reputation

4 Badges

2 years, 342 days
Contact Dkunb

MaplePrimes Activity


These are questions asked by Dkunb

ode system is singular at the initial po...

Dkunb 65 Maple 2021
ode dsolve numeric
April 13 2023
2 1

Could you help me how to deal with this problem?

> 

restart;
> 

with(plots):with(plottools):with(DETools):
> 

 

> 

Sys:=diff(T(R),R)=((1-1/R)*(sqrt(1-(alpha/R)^2*(1-1/R))))^(-1),diff(Phi(R),R)=(alpha/R)^2*(sqrt(1-(alpha/R)^2*(1-1/R)))^(-1);

diff(T(R), R) = 1/((1-1/R)*(1-alpha^2*(1-1/R)/R^2)^(1/2)), diff(Phi(R), R) = alpha^2/(R^2*(1-alpha^2*(1-1/R)/R^2)^(1/2))

 (1)
> 

 

> 

inits:=[[T(0)=0.5,Phi(0)=0],[T(0)=0.5,Phi(0)=Pi/4]];

[[T(0) = .5, Phi(0) = 0], [T(0) = .5, Phi(0) = (1/4)*Pi]]

 (2)
> 

K:=dsolve([Sys,op(op(1,inits))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure);

Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point

  
> 

 

Download HW6.m

Root finder for the fourth power polynom...

Dkunb 65 Maple 2021
solve rootfinding
April 09 2023
1 3

Could you help me to solve this problem for the parameter beta?

> 

restart;
> 

e1:= 0.5; e2:=0.2;theta:=5;yeq:=e2;

.5

  

.2

  

5

  

.2

 (1)
> 

f:=(theta*x-1)*(1-x)*(1+beta*x^2)-y;
g:=x/(1+beta*x^2); gs:=unapply(g,x);

(5*x-1)*(1-x)*(beta*x^2+1)-y

  

x/(beta*x^2+1)

  

proc (x) options operator, arrow; x/(beta*x^2+1) end proc

 (2)
> 

fs:=subs(y=yeq,f);

(5*x-1)*(1-x)*(beta*x^2+1)-.2

 (3)
> 

assumptions:=x>1/theta, x<1,beta>0,beta<1,gs(x)>e1;
solve(fs=0,x,useassumptions) assuming assumptions;

1/5 < x, x < 1, 0 < beta, beta < 1, .5 < x/(beta*x^2+1)

  

Warning, solve may not respect assumed property 'real' on 'x/(beta*x^2+1)'.

  

Error, (in type/realcons) too many levels of recursion

  
> 

gs2:=subs(beta=0.6,gs(x));

x/(.6*x^2+1)

 (4)
> 

sol:= solve(subs(beta=0.6,fs=0),x,useassumptions) assuming x>1/theta, x<1;

.2514854589, .9665623271

 (5)
> 

subs(x=sol[1],gs2);

.2422912423

 (6)
> 

 

Download Rootsfind.mw

animate for the multiple parameters...

Dkunb 65 Maple
textplot plotting animate
April 09 2023
1 1

How can I show all the parameters in title ?

> 

restart;
> 

with(DEtools):with(plots):with(plottools):
> 

 

> 

 

> 

sigma1:=e1*alpha: sigma2:=e2*delta:
> 

g:=x/(1+beta*x^2);
f:=(theta*x-1)*(1-x)*(1+beta*x^2)-y;
h:=alpha*g-z-sigma1;
j:=delta*y-sigma2;

x/(beta*x^2+1)

  

(theta*x-1)*(1-x)*(beta*x^2+1)-y

  

alpha*x/(beta*x^2+1)-e1*alpha-z

  

-delta*e2+delta*y

 (1)
> 

 

> 

 

> 

p0:=theta->plot([1/theta,y,y=0..1],linestyle=dash,color= green):
p1:=e1->plot([x,e1,x=0....1.5],color=blue):
q0:=animate(p0,[theta],theta=2...10):
q1:=animate(p1,[e1],e1=0.1..1):
q2:=plot([1,y,y=0..1],linestyle=dash,color= green):
p3:=beta->plot([x,x/(1+beta*x^2),x=0..1.5],color=magenta);
q3:=animate(p3,[beta],beta=0..1.5):

display([q0,q1,q2,q3],view=[0..1.5,0..1]);

proc (beta) options operator, arrow; plot([x, x/(1+beta*x^2), x = 0 .. 1.5], color = magenta) end proc

  

 
> 

 

Download animate_plots.mw

CharacteristicPolynimial...

Dkunb 65 Maple
spelling
April 07 2023
1 1

I could not figure how to fix my code. Could you help me?

Thanks in advance,

LinearAnalysis.mw

> 

restart;
> 

:with(DynamicSystems):with(DEtools):with(plots):with(plottools):with(DETools):
> 

wiht(LinearAlgebra): with(VectorCalculus):
> 

 

> 

 

> 

e1:=0.5: e2:=0.2:
> 

sigma1:=e1*alpha: sigma2:=e2*delta:
> 

g:=x/(1+beta*x^2);
f:=(theta*x-1)*(1-x)*(1+beta*x^2)-y;
h:=alpha*g-z-sigma1;
j:=delta*y-sigma2;
F:=g*f;
H:=y*h;
G:=z*j;

x/(1+.4*x^2)

  

(5*x-1)*(1-x)*(1+.4*x^2)-y

  

x/(1+.4*x^2)-z-.5

  

y-.2

  

x*((5*x-1)*(1-x)*(1+.4*x^2)-y)/(1+.4*x^2)

  

y*(x/(1+.4*x^2)-z-.5)

  

z*(y-.2)

 (1)
> 

 

> 

theta:=5;alpha:=1;delta:=1;beta:=0.4;
Ffunc:=unapply(F,x,y); Hfunc:=unapply(H,x,y,z);Gfunc:=unapply(G,y,z);
J1:=Jacobian([Ffunc(x,y),Hfunc(x,y,z),Gfunc(y,z)],[x,y,z]);
C:= CharacteristicPolynomial(J1,lambda);
R:= RouthTable(C,lambda);

theta := 5

  

alpha := 1

  

delta := 1

  

beta := .4

  

Ffunc := proc (x, y) options operator, arrow; x*((5*x-1)*(1-x)*(1+.4*x^2)-y)/(1+.4*x^2) end proc

  

Hfunc := proc (x, y, z) options operator, arrow; y*(x/(1+.4*x^2)-z-.5) end proc

  

Gfunc := proc (y, z) options operator, arrow; z*(y-.2) end proc

  

Matrix(3, 3, {(1, 1) = ((5*x-1)*(1-x)*(1+.4*x^2)-y)/(1+.4*x^2)-.8*x^2*((5*x-1)*(1-x)*(1+.4*x^2)-y)/(1+.4*x^2)^2+x*((5-5*x)*(1+.4*x^2)-(5*x-1)*(1+.4*x^2)+(4.0*x-.8)*(1-x)*x)/(1+.4*x^2), (1, 2) = -x/(1+.4*x^2), (1, 3) = 0, (2, 1) = y*(1/(1+.4*x^2)-.8*x^2/(1+.4*x^2)^2), (2, 2) = x/(1+.4*x^2)-z-.5, (2, 3) = -y, (3, 1) = 0, (3, 2) = z, (3, 3) = y-.2})

  

Error, (in Verify:-CommonExports) system object is not a module

  

Matrix(%id = 36893491070243237028)

 (2)
> 

E0:= [0,0,0];
JE0:=subs(x=E0[1],y=E0[2],z=E0[3],J1);
CE0:=CharacteristicPolynomial(JE0,lambda);
RE0:= RouthTable(CE0,lambda);
RouthTable(CE0,lambda,'stablecondition'=true);

E0 := [0, 0, 0]

  

Matrix(3, 3, {(1, 1) = -1., (1, 2) = -0., (1, 3) = 0, (2, 1) = 0., (2, 2) = -.5, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -.2})

  

Error, (in Verify:-CommonExports) system object is not a module

  

Matrix(%id = 36893491070243233172)

  

true

 (3)
> 

E1:= [1,0,0];
JE1:=subs(x=E1[1],y=E1[2],z=E1[3],J1);
CE1:=CharacteristicPolynomial(JE1,lambda);
RE1:= RouthTable(CE1,lambda);
RouthTable(CE1,lambda,'stablecondition'=true);

E1 := [1, 0, 0]

  

Matrix(3, 3, {(1, 1) = -4.000000000, (1, 2) = -.7142857143, (1, 3) = 0, (2, 1) = 0., (2, 2) = .2142857143, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -.2})

  

Error, (in Verify:-CommonExports) system object is not a module

  

Matrix(%id = 36893491070243221500)

  

true

 (4)
> 

E2:= [1/theta,0,0];
JE2:=subs(x=E2[1],y=E2[2],z=E2[3],J1);
CE2:=CharacteristicPolynomial(JE2,lambda);
RE2:= RouthTable(CE2,lambda);
RouthTable(CE2,lambda,'stablecondition'=true);

E2 := [1/5, 0, 0]

  

Matrix(3, 3, {(1, 1) = .8000000000, (1, 2) = -.1968503937, (1, 3) = 0, (2, 1) = 0., (2, 2) = -.3031496063, (2, 3) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = -.2})

  

Error, (in Verify:-CommonExports) system object is not a module

  

Matrix(%id = 36893491070243218004)

  

true

 (5)
> 

 

Download LinearAnalysis.mw

 

Check positivity...

Dkunb 65 Maple 2021
assumptions assuming
March 17 2023
1 5

Why did the last line is(ysol2[1]>0) give false?

Maple does not recognize the assumption?

Thanks in advance.
 

> 

restart;
> 

 

> 

 

> 

interface(showassumed=0);

0

 (1)
> 

assume(theta>1,alpha>0,sigma1>0,beta>0,sigma2>0,delta>0,x>0,y>0,z>0);
> 

 

> 

f:=((theta*x-1)*(1-x)-y);
g:=y/(1+beta*y^2);
h:=(alpha*x-sigma1)*(1+beta*y^2);
j:=(delta*g-sigma2);
dxdt:=x*f;
dydt:=g*(h-z);
dzdt:=z*j;

(theta*x-1)*(1-x)-y

  

y/(beta*y^2+1)

  

(alpha*x-sigma1)*(beta*y^2+1)

  

delta*y/(beta*y^2+1)-sigma2

  

x*((theta*x-1)*(1-x)-y)

  

y*((alpha*x-sigma1)*(beta*y^2+1)-z)/(beta*y^2+1)

  

z*(delta*y/(beta*y^2+1)-sigma2)

 (2)

case0; x = 0, y = 0, z = 0

> 

E0:=<0,0,0>;

Vector[column](%id = 36893490685786064044)

 (3)

"case1: z=0 , h(x,y)=0, f(x,y)=0 : From y=(theta*x-1)*(1-x)=y -> 1/(theta)<x and x<1"

> 

xsol1:=solve(h=0,x) assuming x>1/theta and x<1;
ysol1:=solve(subs(x=xsol1,f)=0,y)assuming xsol1>1/theta and xsol1<1;
E1:=<xsol1,ysol1,0>;

xsol1 := sigma1/alpha

  

ysol1 := -(-sigma1*theta+alpha)*(alpha-sigma1)/alpha^2

  

Vector[column](%id = 36893490685786055012)

 (4)

case2; g(y) = 0, h(x, y) = 0

> 

ysol:=[solve(j=0,y)];
ysol2:=simplify(subs(sigma2=delta*eta,ysol)) assuming (1-4*beta*eta^2>0 and eta>0);

Warning, solve may be ignoring assumptions on the input variables.

  

[(1/2)*(delta+(-4*beta*sigma2^2+delta^2)^(1/2))/(beta*sigma2), -(1/2)*(-delta+(-4*beta*sigma2^2+delta^2)^(1/2))/(beta*sigma2)]

  

[(1/2)*(1+(-4*beta*eta^2+1)^(1/2))/(beta*eta), (1/2)*(1-(-4*beta*eta^2+1)^(1/2))/(beta*eta)]

 (5)
> 

is(ysol2[1]>0)

false

 (6)
> 

 

(1/2)*(1+(-4*beta*eta^2+1)^(1/2))/(beta*eta)

 (7)
> 

 

> 

 


 

Download Project_v1.mw

1 2 3 4 5 Page 3 of 5