Items tagged with ode ode Tagged Items Feed

Hello,

I have two regimes. Each regime is characterized by system of state differential equations (diff(S(t), t), diff(K(t), t)) and co-state differential equations (diff(psi[S](t), t), diff(psi[Iota](t), t)) as follows: 

(1)

diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))-upsilon

diff(psi[S](t), t) = eta*K(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2), diff(psi[Iota](t), t) = eta*S(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2)

(2)

diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))

diff(psi[S](t), t) = eta*K(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2), diff(psi[Iota](t), t) = eta*S(t)^2*(psi[S](t)-psi[Iota](t))/(w*N*(S(t)+K(t))^2)

The first regime is employed from 0 to t1 (where t1 is unknown) and then regime 2, from t1 to T. I know the initial values for the state variables of the first system at t=0, that is, S(0)=S0 and K(0)=K0, as well as boundary conditions for the co-state variables for regime 2 at t=T, that is, psi_S(T)=a and psi_I(T)=b. 

I also know that my unknown t1 should satisfy the algebraic equation: -c - psi_S(t1)=0, where psi_S(t1) is the solution of the co-state diff equation of the regime 2 at t=t1. 

My question is: If I assume that the systems have not analytical solution, how can I found unknown t1 numerically? Moreover, asssume that all other parameters, such as T, eta, upsilon and others are given.

Hello,

My code records the values I need, however, I need to implement a modulo of 2*Pi on my result for theta. But this leads to a graph with no plots and I'm not sure how to fix it. Any help is greatly aprreciated! 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 (`θIC`, `ωIC`, n) local a, Mh, Msat, G, e, beta, M, Eqns, ICs, soln; option remember; global `ω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) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric, maxfun = 0, output = array([seq(i, i = 0 .. n*Th, Th)])); plots:-odeplot(soln, [modp(theta(t), 2*Pi), omega(t)/`ωH`], 0 .. n*Th, labels = ["θ(t)","ω(t)/ωH"], axes = boxed, style = plottools:-point, size = [.25, .75]) end proc:

plots:-display(HyperionOrbit(.5, 1.8*`ωH`, 10));

Download Poincare_section_Boyd_plot_fixing_theta.mw

bia Man

Hello,

I have a question about poincare sections. I have this piece of code i need to analyse and I want to use a poincare section in order to so. How could I do it? I am interested in theta and omega. Any help is greatly appreciated! Thank you in advance!

Kind regards,

Gambia Man

with(plots):

a := 1.501*10^9:

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

1876321.326

 

0.3348672330e-5

(1)

HyperionOrbit := proc (`θIC`, `ωIC`) local a, Mh, Msat, G, e, beta, M, Eqns, ICs; global `ωH`, Th, soln; 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) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric); odeplot(soln, [theta(t), omega(t)/`ωH`], 0 .. 5*Th, numpoints = 2000, labels = ["θ(t)","ω(t)/ωH"], axes = boxed, size = [.25, .75]) end proc

``

 

Download New_Poincare_section.mw

Hi everyone, I'm working a problem

Given a dog and a man. At t=0, the position of the man and his dog is (0,0) and (0,h), respectively. Then, the man start moving along Ox with constant speed vm. The dog keep running toward its master with constant speed vd. Describe the movement of the dog?

Therefore, let say the position of the dog is (f1(x),f2(x)), I wrote these:

restart;

with(plots);

MD := proc (h, vm, vd) local x, y, ode, ics, sol;

ode := {(diff(f2(t), t))*(vm*t-f1(t))+(diff(f1(t), t))*f2(t), diff(f2(t), t) = -sqrt(vd^2-(diff(f1(t), t))^2)};

ics := {f1(0) = h, f2(0) = 0};

sol := dsolve(ode union ics);

x := unapply(eval(f1(t), sol), t);

y := unapply(eval(f2(t), sol), t);

plot([x(t), y(t), t = 0 .. 10], scaling = constrained) end proc;

MD(10, 1, 5);

However, "sol" is returned NULL, which means the equation has no solution. I supposed I completed the motion part of the problem correctly. Please help me or point out an another method

 

Hello,

I have an error in my code. I don't know where it came from. Earlier today I loaded this and it worked fine and now an error comes up. Any help is greatly appreciated! Thank you in advance!

Kind regards,

Gambia Man

with(plots):

``

HyperionOrbit := proc (`θIC`, `ωIC`) local a, Mh, Msat, G, e, beta, M, Eqns, Th, ICs, soln, `ωH`; 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); Th := sqrt(4*Pi^2*a^3/(G*(Mh+Msat))); `ωH` := 2*Pi/Th; ICs := xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e))), theta(0) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric); odeplot(soln, [theta(t), omega(t)/`ωH`], 0 .. 5*Th, numpoints = 2000, labels = ["θ(t)","ω(t)/ωH"], axes = boxed) end proc

HyperionOrbit(.5, 1.8*`ωH`)

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

 

Error, (in dsolve/numeric/SC/IVPpoints) unable to store '11.3097335529232/Th' when datatype=float[8]

 

``

``

``

 

Download Poincare_section.mw

Hello,

I have an issue with this code. The issue lies in the warning which is outputed. However, I don't know how to resolve it. Any help would be greatly appreciated! Thank you in advance!

Kind regards,

Gambia Man

with(plots):

a := 1.501*10^9:

Eqns := diff(xH(t), t) = vxH(t), diff(vxH(t), t) = -G*(Mh+Msat)*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(yH(t), t) = vyH(t), diff(vyH, t) = -G*(Mh+Msat)*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*(Mh+Msat)*(1-e)/(a*(1+e))):

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

soln := dsolve({Eqns, ICs}, numeric):

plots:-odeplot(soln, [xH(t)/a, yH(t)/a], 0 .. 10*T, scaling = constrained, labels = ["x/a", "y/a"], numpoints = 2000);

Warning, cannot evaluate the solution further right of 453574.15, probably a singularity

 

 

``

 

Download Hyperion_Orbit.mw

 

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!

 

This is the code hw2_final.mw

This is the warning 

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters )"

 

How to solve it?

Hello,

 I am trying to solve an differential equation in Maple, but I have no response from Maple.Can someone give a look for me? The coefficients a and b cannot be zero.

 Thanks a lot.

 

Download edo_solução_A_U.mw

Hello everybody

I have a system of ODEs that are too large. I am trying to solve this system using Maple. The corrsponding file has been uploaded. I have 2 question. 1- Does anyone know what is the source of error? 2- Is this a high expectation for Maple to solve such a too big system of ODEs. Note: It makes a few minutes to run all the programm.

Thanks in advance

M.mw

Hello

i have an ODE like this:

I sove this ODE with plot order:

with(plots);
odeplot(sol, [x, (3*D1*a+4*D2)*P(x)/((1-q*S(x))*D2)], .5 .. (1/2)*Pi, tickmarks = [[seq((1/10)*i*Pi = (180*i*(1/10))*`°`, i = 1 .. 8)], default]);
my plot work very well. but i need to plot this ODE with five different parameter (q for for instance, q=0.1 & q=0.2 ....) all in one axis. something like this:

plot(sin(x),x=-10..10)hello

I have an ODE plot like this and I want its horizontal axes to be in degree instead of radian, but I don't know how

Hello,

In this procedure I get the error, 

Error, `:=` unexpected

I know what it means I just can't seem to resolve it. Any help would be greatly appreciated! Thank you in advance for looking at this code!

Kind regards,

Gambia Man

HamilMat := proc (K::integer) local ni, mi, nj, mj, N, Hamil, Eigenvec, i, j, res; option remember; global Vij, U, L; N := K^2; ni := Vector(N); mi := Vector(N); nj := Vector[row](N); mj := Vector[row](N); for i to N do for j to K do res := (i+K-j)/K; if type(res, integer) = true then ni[i] := j; nj[i] := j; mi[i] := res; mj[i] := res end if end do end do; Hamil := Matrix(N, shape = symmetric); for i to N do for j from i to N do if i <> j then Hamil(i, j) := Vij(ni[i], mi[i], nj[j], mj[j]) elif i = j then Hamil(i, j) := Vij(ni[i], mi[i], nj[j], mj[j])+(1/2)*(ni[i]^2+mi[i]^2)*Pi^2/L^2 end if end do end do; return Eigenvec := Eigenvectors(Hamil, output = ['values', 'vectors']), Hamil end proc

Error, `:=` unexpected

 

 

Download small_error.mw

 with small step rather than two points?

I had writen my questions into this code hw2_finished.mw And I think there's maybe something wrong with A=[0..130]

I think you can understand what I want by reading this code but if you are confused please let me know.

Hi everyone, I'm a new one to Maple. I've just learnt some basic tools :)

 

This is my task. I tried to record in Maple but I had errors. I don't know why I had problems but I hope you will help me and I will do it.

tkanks

3 4 5 6 7 8 9 Last Page 5 of 29