Items tagged with odeplot odeplot Tagged Items Feed

I have some differential equations that I want to plot on the same axis (as i have below), but would like to plot with a log scale to illustrate the fact that one is simply a logarithmic decay (solid line) and all the others are not.

I had used deplot and display to make the above graph,

but if you want more detail, here is a worksheet with the relevant differntial equations:
MaplePrimesGraphUpload.mw

 

 

 

Hi everyone!

I tried to plot the solution of the following ode, but I only got the message error:

Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

(file attached)

Problem.mw

https://dl.dropboxusercontent.com/u/66502600/Problem.mw

 

Please, help me!

 

Thank you so much!

 



thanks. I played around, and had problems implementing your ideas for one of the systems I'm interested in.I don't see a difference between this and what you had advised me on, but it gets an error.

any idea why?
or how to fix it?

thing1 := diff(B[1](t), t) = piecewise(t <= 500, 0.3e-2-(63/10000)*B[1](t)-(3/500)*B[2](t), -(3/10000)*B[1](t)):
thing2 := diff(B[1](t), t) = piecewise(t <= 500, 0.1e-1-(1/50)*B[1](t)-(13/625)*B[2](t), -(1/1250)*B[2](t)):
sol := dsolve({thing1, thing2, B[1](0) = 0, B[2](0) = 0}, {B[1](t), B[2](t)}, numeric, output = listprocedure); plots:-odeplot(sol, [B[1](t), B[2](t)], t = 450 .. 550);

Error, (in dsolve/numeric/DAE/explicit) unable to obtain the standard form of the DAE system due to the presence of leading dependent variables/derivatives in the piecewise: piecewise(t <= 500, 1/100-(1/50)*B[1](t)-(13/625)*B[2](t), -(1/1250)*B[2](t))-piecewise(t <= 500, 3/1000-(63/10000)*B[1](t)-(3/500)*B[2](t), -(3/10000)*B[1](t))
Error, (in plots/odeplot) curve is not fully specified in terms of the ODE solution, found additional unknowns {B[1](t), B[2](t)}


Hello,

I would like to plot gait diagrams (the lines you can see on the picture belowà from the solutions obtained with a NL oscillator (composed with 8 coupled odes). Here the result that I would like to obtain.

Initial plot:

 

Desired plot

 

 

I would like to obtain 4 lines corresponding to the 4 elliptic trajectories obtained with the NL oscillator. The four lines should be done like this. When the trajectory is above 0, the line should be colored in green. When the trajectory is below 0, the line should be colored in black. 

May you help me to define this kind of graph called gait diagrams from the solution of the NL oscillator ?

Here you can find my maple code:

K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);

for i to 4
do
r[i]:=sqrt((u[i](t))^2+(v[i](t))^2):
omega[i]:=omega[sw]/(1+exp(b*v[i](t)))+omega[st]/(1+exp(-b*v[i](t))):
Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)-omega[i]*v[i](t):
Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*u[i](t)+MatrixVectorMultiply(K,<seq(v[i](t),i=1..4)>)[i]:
EqSys[i]:=[Equ[i],Eqv[i]]:
end do:

paramsCycle:=omega[st]=4*2*Pi,omega[sw]=2*Pi,Au=5,Av=50,b=100;
params:=paramsCycle;

Differential system 
sys:=map(op,eval([seq(EqSys[i],i=1..4)],[params]));
ic:=[u[1](0)=0, v[1](0)=0,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1];
Résolution1
res:=dsolve([sys[],ic[]],numeric):
Initial boundaries
tcalc:=4;
ic2:=[seq(u[i](0)=eval(u[i](t), res(tcalc)),i=1..4),seq(v[i](0)=eval(v[i](t), res(tcalc)),i=1..4)];
Résolution2
res:=dsolve([sys[],ic2[]],numeric):

tmax:= 40:
numpts:=100*tmax:
plots:-odeplot(res,[t,v[1](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[t,v[2](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[t,v[3](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[t,v[4](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);
plots:-odeplot(res,[seq([t,v[i](t)+i*5],i=1..4)],0..tmax,thickness=2,view=[0..5,0..25], numpoints = numpts);

Thanks a lot for your help

Dear all;

Thank you for helping me to solve this  question.

I solve an ode, but I have an error when I would like to plot the solution.

uanble to achieve continuous solution with requested accuracy of 0.1e-5 with maximum 128 point mesh (was able to get 0.14e-5), consider increasing `maxmesh` or using larger `abserr`
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

I try to increase the point mesh or take a large abserr but always I have the same problem.

 

restart:
with(plots):

ode := diff(y(x), x, x) = x*y(x)+sqrt(x);

ics := y(0) = 0, y(1) = 1;
sol:=dsolve({ode,ics}, numeric):
odeplot( sol,[x, y(x)], x=0..1, maxmesh=1000);

hi guys,

 

i have a first order differential equation and i want to plot it with odeplot in polar coordinates .

 

 

thanks in advance1.mw

I am exploring the onset of chaos in the driven damped pendulum. I have set up the non-linear initial value problem that describes the situation in terms of phi(t), the angle of the pendulum from the vertical. I can solve the differential equation numerically for certain initial conditions using dsolve, and then graph that solution using odeplot. I can then explore what happens when I change the initial conditions ever so slightly. Again I can solve this new initial value problem numerically with dsolve, and then graph it with odeplot, Suppose the solution to the first problem is phi1(t), and the solution to the second initial value problem phi2(t). I would like to plot log | phi2(t) - phi1(t) | versus t, that is, the logarithm of the absolute value of the difference of the two solutions versus t. I’m interested in how sensitive the system is to initial conditions.

So, how can I disentangle the two solutions, phi1(t) and phi2(t), from dsolve and then plot the graph I want?


Hello,

I've been trying to collect points at certain intervals in my plot, however, the sample option doesn't seem to be working. Can anyone lend me a hand? Any help is greatly appreciated! Thank you in advance!

Kind regards,

Gam

with(plots):

a := 1.501*10^9:

Th := sqrt(4*Pi^2*a^3/(G*(Mh+Msat))):

HyperionOrbit := proc (`&theta;IC`, `&omega;IC`, n) local a, Mh, Msat, G, e, beta, M, Eqns, ICs, soln; global `&omega;H`, Th; a := 1.501*10^9; Mh := 5.5855*10^18; Msat := 5.6832*10^26; G := 6.67259/10^11; e := .232; beta := .89; M := Mh+Msat; Eqns := diff(theta(t), t) = omega(t), diff(omega(t), t) = -G*Msat*beta^2*(xH(t)*sin(theta(t))-yH(t)*cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t)))/(xH(t)^2+yH(t)^2)^2.5, diff(xH(t), t) = vxH(t), diff(vxH(t), t) = -G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(yH(t), t) = vyH(t), diff(vyH(t), t) = -G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2); ICs := xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e))), theta(0) = `&theta;IC`, omega(0) = `&omega;IC`; soln := dsolve({Eqns, ICs}, numeric, maxfun = 0); plots:-odeplot(soln, [theta(t), omega(t)/`&omega;H`], 0 .. n*Th, numpoints = 2000, labels = ["&theta;(t)","&omega;(t)/&omega;H"], axes = boxed, style = point, sample = [seq(i, i = 0 .. n*Th, Th)]) end proc:

HyperionOrbit(.5, 1.8*`&omega;H`, 10)

Error, (in plot/options2d) unexpected option: sample = [0, 1876321.326, 3752642.652, 5628963.978, 7505285.304, 9381606.630, 11257927.96, 13134249.28, 15010570.61, 16886891.93, 18763213.26]

 

``

``

 

Download Poincare_section.mw

bia Man

Suppose we have a differential equation that it's order is 2. For example

diff(x(t), t, t)+(1000*(x(t)^2-1))*(diff(x(t), t))+x(t) = 0,

and we want to solve it numerically. When we use some initial values and dsolve it, the maple solves it very quickly. We can plot each new equations of x(t) without any problem respect of "t" or else.

But when we want to show the variation of "diff(x(t), t, t)" respect of "t" (when our derivation order is equal or upper than our differential equation order) our answer is 0 however dx(t)/dt is not 0 or constant!

What's happen here?

I tested it with all kind of methods as "rkf45","rosenbrock","lsode" etc.  but non of them showed me a correct answer.

How can I solve it correctly?
 
  
Please help me!

 

Hello,

I'm writing a code and I seem to have an issue when trying to implement a procedure. Here is the code:

 

with(plots):

Z := 75; A := 189; k := 14.6; Rm := 8*R; r0 := 10^(-8)*R; c := 137.036; ms := 105.66/(.51100)

fmtoau := 10^(-15)/(0.529177e-1*10^(-9)):

R := 1.1*fmtoau*A^(1.0/(3.0)):

f := proc (x) options operator, arrow; 1/(1+exp(k*(x-1))) end proc:

n0 := 3*Z*k^2/(4*Pi*(Pi^2+k^2)*R^3):

n := proc (r) options operator, arrow; 4*Pi*n0*f(r/R) end proc:

int(r^2*n(r), r = 0 .. Rm);

73.00000007

(1)

plot(n/n0, 0 .. 2*R);

 

v1 := unapply(int(x^2*f(x), x), x):

Vfermi := proc (r) options operator, arrow; -4*Pi*n0*R^2*(R*(v1(r/R)-v10)/r+v2Rm-v2(r/R)) end proc:

Vuniform := proc (r) options operator, arrow; piecewise(r < R, -Z*(3/2-(1/2)*r^2/R^2)/R, -Z/r) end proc:

plot([Vuniform(r), Vfermi(r), -Z/r], r = r0 .. 2*R, V = 1.2*Vfermi(r0) .. 0, legend = ["uniform charge", "Fermi distribution", "point charge"]);

 

``

plotsol1s := proc (E, K, r0, S, col) local Eqns, ICs, fnl, gnl, r, soln; global ms, c; Eqns := diff(fnl(r), r) = gnl(r)*[E+2*ms*c^2-Vuniform(r)]/c-(K+1)*fnl(r)/r, diff(gnl(r), r) = -fnl(r)*[E-Vuniform(r)]/c-(1-K)*fnl(r)/r; ICs := fnl(r0) = 1, gnl(r0) = 0; soln := dsolve({Eqns, ICs}, numeric); plots:-odeplot(soln, [r, fnl(r)], r0 .. S, color = col) end proc:

plotsol1s(-3*10^5, -1, 10^(-10), Rm, red)

Error, (in f) unable to store '[HFloat(0.0)]' when datatype=float[8]

 

``

``


Any help would be greatly appreciated.

 

Gambia Man

 

sys := [(diff(c(t), t))*(diff(a(t), t))*(diff(b(t), t))+(diff(c(t), t))*t^2, (diff(c(t), t))*(diff(a(t), t))*t+(diff(c(t), t))*t*(diff(b(t), t)), (diff(c(t), t))*(diff(a(t), t))*(diff(b(t), t))+(diff(c(t), t))*(diff(a(t), t))*t+(diff(c(t), t))*t^2]
DEplot(sys, [a(t), b(t), c(t)], t = 0 .. 2, a = -15 .. 15, b = -15 .. 15, c = -15 .. 15, color = magnitude, title = `Stable Limit Cycles`, arrows = curve, dirfield = 800, axes = none);
odeplot(sys, [t, a(t), b(t), c(t)], -4 .. 4, color = orange);

restart;
odes:=diff(x1(t),t)*diff(x2(t),t$2)*sin(x1(t)*x2(t))+5*diff(x1(t),t$2)*diff(x2(t),t)*cos(x1(t)^2)+t^2*x1(t)*x2(t)^2=exp(-x2(t)^2),diff(x1(t),t$2)*x2(t)+diff(x2(t),t$2)*diff(x1(t),t)*sin(x1(t)^2)+cos(diff(x2(t),t$2)*x2(t))=sin(t);
ics:=x1(0)=1,D(x1)(0)=1,x2(0)=2,D(x2)(0)=2;
subs(diff(x1(t),t$2)=yp2,diff(x2(t),t$2)=yp4,diff(x1(t),t)=Y[2],diff(x2(t),t)=Y[4],x1(t)=Y[1],x2(t)=Y[3],{odes});
p:=proc(N,t,Y,YP)
local eqs,yp2,yp4;
YP[1]=Y[2];
YP[3]=Y[4];
eqs:=[yp2*Y[3]+yp4*Y[2]*sin(Y[1]^2)+cos(yp4*Y[3]) = sin(t), Y[2]*yp4*sin(Y[1]*Y[3])+5*yp2*Y[4]*cos(Y[1]^2)+t^2*Y[1]*Y[3]^2 = exp(-Y[3]^2)];
YP[2],YP[4]:=op(subs(fsolve(eqs,{yp2=1,yp4=2}),[yp2,yp4]));
end proc:
res:=dsolve(numeric,procedure=p,initial=array([1,1,2,2]),number=4,procvars=[x1(t),diff(x1(t),t),x2(t),diff(x2(t),t)],start=0,maxfun=0);
plots:-odeplot(res,[t,x1(t)],0..5,gridlines=true);

It's a long time to wait for the odeplot.

Any advice is appreciated.

Hi all,

I would like to know how could I plot an ODE for 3 values of a parameter in the same figure.

sol:=dsolve(dsys,numeric,parameters=[t]);

p:=proc(t)sol(parameters=[t]); plots:-odeplot(sol,[x,f(x)],-1..0,_rest) end proc;

p(-1); p(-0.5); p(-0.1);

Here, I have 3 figures one for each value. I used with(plots): display(p(-1); p(-0.5); p(-0.1)). But doesn't work.

Thanks for your help.

 

Hello,

After solving ode I am looking only for the values >=1.5. For example at t=1, y(t)=3.8940.

How can I extract the values >= 1.5 from the solution to use it as data (t,y(t)) and save it ?

restart;
with(DEtools); with(plots);
eqn := diff(y(t), t) = -.25*y(t);

 init := y(0) = 5;

sol := dsolve({eqn, init}, {y(t)}, numeric, output = array([seq(i, i = 0 .. 50)]));
p[1] := plot(1.5, t = 0 .. 50, colour = black);

p[2] := odeplot(sol, [t, y(t)], t = 0 .. .50, colour = red);

display(p[1], p[2]);

 

Thanks

Hello,

The idea: parameter "a" will have a new random value each 10 days.

The way I did it is working but it can get very long especially if I do it for a system of equations and for long time more than a year (365 days).

The code:

with(DEtools); with(plots);
n := 5;

for i to n do Ra[i] := RandomTools:-Generate(distribution(Uniform(0.1e-1, .5))); a[[i]] := Ra[i] end do;

b := 0.1e-2;

T := 10;

 eq := diff(L(t), t) = a*L(t)-b;

init[1] := L(0) = 100;
 sol[1] := dsolve({init[1], subs(a = a[[1]], eq)}, L(t), range = 0 .. T, numeric);


init[2] := L(T) = rhs(sol[1](T)[2]);

sol[2] := dsolve({init[2], subs(a = a[[2]], eq)}, L(t), range = T .. 2*T, numeric);

 

init[3] := L(2*T) = rhs(sol[2](2*T)[2]);
sol[3] := dsolve({init[3], subs(a = a[[3]], eq)}, L(t), range = 2*T .. 3*T, numeric);

p[1] := odeplot(sol[1], [t, L(t)], t = 0 .. T);

p[2] := odeplot(sol[2], [t, L(t)], t = T .. 2*T);

p[3] := odeplot(sol[3], [t, L(t)], t = 2*T .. 3*T);

p := display([p[1], p[2], p[3]]);
display(p);

 

Thank you

1 2 3 Page 1 of 3