Items tagged with ode ode Tagged Items Feed

Hi everyone,

I have been  trying to solve a coupled system of 2 differencial equations, 1 PDEs with 1 ODE.

The code is below.

> restart;

> with(linalg):with(plots):
> PDE:=[diff(x(z,t),t)=(a/Pe)*diff(x(z,t),z$2)-a*diff(x(z,t),z)-(1-epsilon)/epsilon*3*Bi*(x(z,t)-1)/(1-Bi*(1-1/ksi(z,t)))]:
> a=2:Pe=3:Bi=5:epsilon=0.85:

> ODE := [diff(ksi(z,t),t) = (b*Bi*x(z,t)-1)/ksi(z,t)^2/(1-Bi*(1-1/ksi(z,t)))]:
> IC1:=c(z,0)=0:
> IC2:=ksi(z,0)=1:
> bc2:=diff(x(z,t),z):
> bc1:=x(z,t)-1/pe*diff(x(z,t),z):
> N:=10:
> L:=1:
> dyduf:=1/2*(-x[2](t)-3*x[0](t)+4*x[1](t))/h:
> dydub:=1/2*(-x[N-1](t)+3*x[N+1](t)+4*x[N](t))/h:
> dydu:=1/2/h*(x[m+1](t)-x[m-1](t)):
> d2ydu2:=1/h^2*(x[m-1](t)-2*x[m](t)+x[m+1](t)):
> bc1:=subs(diff(x(z,t),z)=dyduf,x(z,t)=x[0](t),z=1,bc1):
> bc2:=subs(x(z,t)-1/pe*diff(x(z,t),z)=dydub,x(z,t)=x[N+1](t),t=0,bc2):
> eq[0]:=bc1:
> eq[N+1]:=bc2:
> eq[m]:=subs(diff(x(z,t),z$2)=d2ydu2,diff(x(z,t),z)=dydu,diff(x(z,t),t)=dydu,x(z,t)=x[m](t),z=m*h,PDE):
> for i from 1 to N do eq[i]:=subs (m=i,eq[m]);od:
> x[0](t):=(solve(eq[0],x[0](t)));

> x[N+1](t):=solve(eq[N+1],x[N+1](t));

> for i from 1 to N do eq[i]:=eval(eq[i]);od:
> eqs:=[seq((eq[i]),i=1..N)]:
> Y:=[seq(x[i](t),i=1..N)];

> A:=genmatrix(eqs,Y,'B1'):
Error, (in linalg:-genmatrix) equations are not linear

> evalm(B1);

[B1[1], B1[2], B1[3], B1[4], B1[5], B1[6], B1[7], B1[8], B1[9],


> B:=matrix(N,1):for i to N do B[i,1]:=B1[i]:od:evalm(B);

> h:=eval(L/(N+1));


> A:=map(eval,A);


> if N > 10 then A:=map(evalf,A);end:
> evalm(A);


> mat:=exponential(A,t);
Error, (in linalg:-matfunc) input must be a matrix

> mat:=map(evalf,mat):
> Y0:=matrix(N,1):for i from 1 to N do
> Y0[i,1]:=evalf(subs(x=i*h,rhs(IC)));od:evalm(Y0);
Error, invalid input: rhs received IC, which is not valid for its 1st argument, expr


> s1:=evalm(Y0+inverse(A)&*b):
Error, (in linalg:-inverse) expecting a matrix

> Y:=evalm(mat&*s1-inverse(A)&*b):
Error, (in linalg:-inverse) expecting a matrix

> Y:=map(simplify,Y):
> Digits:=5;

Digits := 5

> for i from 1 to N do x[i](t):=evalf((Y[i,1]));od:
Error, invalid subscript selector

> for i from 0 to N+1 do x[i](t):=eval(x[i](t));od;


> setcolors(["Red", "Blue", "LimeGreen", "Goldenrod", "maroon",
> "DarkTurquoise", "coral", "aquamarine", "magenta", "khaki", "sienna",
> "orange", "yellow", "gray"]);

["Red", "LimeGreen", "Goldenrod", "Blue", "MediumOrchid",


> pp:=plot([seq(x[i](t),i=0..N+1)],t=0..10,thickness=4);pt:=textplot([[0.28,0.15,typeset("Follow the arrow: ",x[0],"(t), ..., ",
> x[N+1],"(t).")]]):
Warning, unable to evaluate the functions to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct


VIEW(0. .. 10., DEFAULT))

> display([pp,pt,arw],title="Figure 5.13",axes=boxed,labels=[t,"x"]);
Error, (in plots:-display) expecting plot structures but received: [arw]


In advance, thanks for the time of reading it!




I want to solve two odes with their boundary condition. I wrote the code below:




sys_ode:=  eq2=0,eq3=0;
bcs := phi(0)=0,phi(h)=1,T(0)=0,T(h)=1;
sol:=dsolve([sys_ode, ics]);

however, this code doesnt get my desired results (the results are complex!). but when I (with hand) integrate Eq3 twice and substitute boundary conditions and replace in Eq2 the answer is easy and straightforward.

How can I change the following algorithm to get my results?

Thanks for your attention in advance


I am trying to solve the rate equations for A<=>B<=>C with k1,k2,k3,k4 as the rate constants.  I put in the equations into Maple but can't seem to get it to solve them.  The initial conditions for A(0) = S, B(0)=C(0)=0


Here is what I plugged into the program:


Parameters(S, k1, k2, k3, k4);

sys_ode := d*A(t)/dt = -k1*A(t)+k2*B(t), d*B(t)/dt = k1*A(t)+(-k2-k3)*B(t)+k4*C(t), d*C(t)/dt = k3*B(t)-k4*C(t);

d A(t)
------ = -k1 A(t) + k2 B(t),

d B(t)
------ = k1 A(t) + (-k2 - k3) B(t) + k4 C(t),

d C(t)
------ = k3 B(t) - k4 C(t)

ics := A(0) = S, B(0) = 0, C(0) = 0;
A(0) = S, B(0) = 0, C(0) = 0

dsolve({ics, sys_ode});
Error, (in dsolve) found functions with same name but depending on different arguments in the given DE system: . It is required an indication of the dependent variables

I encounter "insufficient  initial/boundary value" error message,  do know how to proceed from there, search with "insufficient initial value" gets no result. Any help will be appreciated.



> restart; alias(r = r(t), f = f(t)); with(plots);
r, f
> DE := diff(r, t) = 2*r+alpha*r*f, diff(f, t) = -f+alpha*r*f;
d d
--- r = 2 r + alpha r f, --- f = -f + alpha r f
dt dt
> params := alpha = .3;
alpha = 0.3
> initv := r(0) = 101, f(0) = 2;
> dvars := [r, t];
> chaodisplay := proc (chartname) EQ := [op(subs(params, [DE])), initv]; EQ1 := dsolve(EQ, numeric); odeplot(EQ1, dvars, t = 0 .. 300, axes = frame, numpoints = 50000, color = green, orientation = [-30, 100], title = chartname) end proc;
Warning, `EQ` is implicitly declared local to procedure `chaodisplay`
Warning, `EQ1` is implicitly declared local to procedure `chaodisplay`
> chaodisplay("Rabbit and Fox");
Error, (in dsolve/numeric/type_check) insufficient initial/boundary value information for procedure defined problem

Dear experts;

How can I solve this problem with maple?



ode:=diff(x[1](t), t)=(lambda-d*x[1](t)-(1-beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1])*beta*x[1](t)*x[3](t)),
 diff(x[2](t), t) =((1-beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1])*beta*x[1](t)*x[3](t)-delta*x[2](t)),
 diff(x[3](t), t) =((1+psi[3](t)*p*x[2](t)/A[2])*p*x[2](t)-c*x[3](t)),diff(psi[1](t), t) =-1+1/A[1]*beta^2*x[1](t)*(x[3](t))^2*(psi[1](t)-psi[2](t))^2-psi[1](t)*(-d+beta^2*(x[3](t))^2*(psi[1](t)-psi[2](t))/A[1]*x[1](t)-(1-beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1])*beta*x[3](t))-psi[2](t)*(-beta^2*(x[3](t))^2*(psi[1](t)-psi[2](t))/A[1]*x[1](t)+(1-beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1])*beta*x[3](t)),
 diff(psi[2](t), t) =1/A[2]*psi[3](t)^2*p^2*x[2](t)+psi[2](t)*delta-psi[3](t)*(psi[3](t)*p^2/A[2]*x[2](t)+(1+psi[3](t)*p*x[2](t)/A[2])*p),
 diff(psi[3](t), t) = 1/A[1]*beta^2*(x[1](t))^2*x[3](t)*(psi[1](t)-psi[2](t))^2-psi[1](t)*(beta^2*(x[1](t))^2*(psi[1](t)-psi[2](t))/A[1]*x[3](t)-(1-beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1])*beta*x[1](t))-psi[2](t)*(-beta^2*(x[1](t))^2*(psi[1](t)-psi[2](t))/A[1]*x[3](t)+(1-beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1])*beta*x[1](t))+psi[3](t)*c;

ics := x[1](0)=5.5556*10^7, x[2](0)=1.1111*10^7,x[3](0)=6.3096*10^9,psi[1](100)=0,psi[2](100)=0,psi[3](100)=0;

dsolve([ode, ics],numeric);?????????????????????????

Please help me


Minimize doesn't work with dsolve porcedure?

tr := proc (x, y)::integer; tr := x+y; result := x^(2+y) end proc

Warning, `result` is implicitly declared local to procedure `tr`


tr(5, 5)



with(Optimization); Minimize(tr(x, y), x = 0 .. 1000, y = 1 .. 1, initialpoint = {x = 25, y = 1})

[0.117556065072605623e-15, [x = HFloat(4.898709434833346e-6), y = HFloat(1.0)]]


xxx := 97.39391293; yyy := -1.588898710



xx := 100;



trool := proc (leng, alpha)::integer; global psi, zx, zy, xx, yy, xxx, yyy, sa, ca, ps, Vx, Vy, vx, vy, ode, ics, XX, YY, trool, G, str, start, ds; sa := evalf(sin(alpha)); ca := evalf(cos(alpha)); ps := evalf(evalc(Im(evalc(str*(x+I*y)-((1/2)*I)*G*ln(x+I*y-start)/Pi)))); psi := ps; xxx := evalf(xx+leng*ca); yyy := evalf(yy+leng*sa); Vx := diff(psi, y); Vy := -(diff(psi, x)); vx := Re(evalf(subs(x = xxx, y = yyy, subs(vvx = Vx, vvx)))); vy := Re(evalf(subs(x = xxx, y = yyy, subs(vvy = Vy, vvy)))); proc (X) options operator, arrow; X(t) end proc; proc (Y) options operator, arrow; Y(t) end proc; zx := proc (t) options operator, arrow; evalf(subs(x = X(t), y = Y(t), subs(vvx = Vx, vvx))) end proc; zy := proc (t) options operator, arrow; evalf(subs(x = X(t), y = Y(t), subs(vvy = Vy, vvy))) end proc; ode := diff(X(t), t) = zx(t), diff(Y(t), t) = zy(t); ics := X(0) = xxx, Y(0) = yyy; ds := dsolve([ode, ics], type = numeric, [X(t), Y(t)], method = rkf45, maxfun = 0, output = listprocedure, abserr = 0.1e-3, relerr = 0.1e-3, minstep = 0.1e-1); XX := rhs(ds[2]); YY := rhs(ds[3]); trool := XX(0.1e-3) end proc:

with(Optimization); Minimize(trool(alpha, leng), assume = nonnegative, alpha = 0 .. 2*Pi, leng = .2 .. 2, iterationlimit = 1000, initialpoint = {alpha = 1, leng = 1})

Error, (in XX) parameter 'alpha' must be assigned a numeric value before obtaining a solution


alpha = 0 .. 2*Pi, leng = .2 .. 2, output = solutionmodule

alpha := 1; leng := 1; XX(10)




trool(1, 11)









with trool procedure minimize dosent work .... and its make me realy sad, couse i need to optimize alpha and leng in other (big one) porcedure with same dsolve.

get this errors:
"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 XX) parameter 'alpha' must be assigned a numeric value before obtaining a solution"

I have following expression


which is 1 solution of the ODE

ode2 := -(diff(y(t), t, t))+(4-12/(1+s*cosh(2*t))+(8*(-s^2+1))/(1+s*cosh(2*t))^2)*y(t) = 0

Now I wanted to construct 2 linear independent solutions via:



and calculate the Wronskian:

with(LinearAlgebra); with(VectorCalculus)

Determinant(Wronskian([f(t_b-t), f(t-t_a)], t))

Since I know these functions are solutions of the second order ODE which does not contain any first order derivative the Wronskian should be a constant. Unfortunately Maple has a hard time to simplify it since the epxression is a little big. Is it my fault or has anyone an idea what to do?

I have a problem in excuting this differential equation in maple it takes a long time but yet no result.

> restart;

> Delta:= epsilon[2]-epsilon[1];

> epsilon[y] := epsilon[2]-(1/4)*Delta*(1-tanh(a*y))^2;

 > z:= tanh(a*y) ;

 > ODE[4]:= diff(Y(y),y,y)- ( a/2* Delta *(1-z)*(z^2-1))/(epsilon[2]- Delta*(1-z)/4)* diff(Y(y),y)-( beta^2+ mu[0]*epsilon[y]*omega^2)*Y(y) = 0;

> dsolve(ODE[4],Y(y));

does this always occur or i do have problem with my version of maple 15, 7 and 16.

Thank you, looking forward for your answers.


I am trying to solve a nonlinear second order ODE with a parameter to be determined. I can set up most of the problem but I am having trouble trying to tell the computer the following boundary condition,


diff(f(x),x) = -K on f(x)=0.


(K will be inputted and is not to be solved for) As i said before, the other boundary conditions are fine and the numerical solution works if i use different boundary conditions. 


For other boundary conditions (for example df/dx = 0 at x=0) I write in the form 



I hope this makes sense and someone has a solution. Thanks in advance.



Given the following system of first order ODE,

dx/dt=0.2x(1-0.5 x)-(1.5 xy)/(1+0.116 x),

dy/dt=(1.3 xy)/(1+0.1x)-0.8y.


Draw a DEplot (for t from 0 to 50) and indicate the particular

trajectory with the initial conditions x(0)=1,y(0)=2. If I

switched to forward Euler method,what would the DE plot look

like then? Is it possible to make the plot made by the

forward Euler method look close to the one which used the

default method?

As it says in the title, I would like to solve the following ODE numerically using forward Euler method, without using the Student Package.



with initial condition y(0)=1. I want to solve it for up to t=1, and then plot both the solution by Euler's method and the solution by "dsolve" on the same graph so I can compare them.


Also, can I make a separate DEplot with t extending to 5?


Thanks in advance.

Following previous question at

and also

I wrote the following code






eq2:=diff(T(eta),eta,eta)+1/(k(eta)/k1[w])*(2/(1-zet^2)*rho(eta)*c(eta)*u(eta)/(p2*10000)+( (a[k1]+2*b[k1]*phi(eta))/(1+a[k1]*phi1[w]+b[k1]*phi1[w]^2)*diff(phi(eta),eta)-k(eta)/k1[w]/(1-eta)*diff(T(eta),eta) )):
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):



p:=proc(pp2) global res,F0,F1,F2:
if not type([pp2],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({eq1=0,subs(p2=pp2,eq2)=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=phi0,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure):
evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2*10000:
end proc:




phb:=evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet))) / evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet))) :



as you can see at the second line of the code, the value of phi0:=0.00789. however, I want to modify the code in a way that phi0 is calculated with the following addition constraint

evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet))) / evalf(2/(1-zet^2)*(Int((1-eta)*F0(eta),eta=0..1-zet)))-0.02=0

I would be most grateful if you could help me in this problem.

Thanks for your attention in advance


Here is a screen of the original question

 My advice to the questioner is to visit a psychiatrist ASAP.

Markiyan Hirnyk


Thank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank youThank you

So i have to plot 3 ODEs. How would I go about plotting dx/dt=10(y-x) AND dy/dt=x(28-x)-y AND dz/dt= xy-8z/3. Any help at all would be appreciated.

The question is to solve

x^2*(diff(y(x), x, x))+x^3*(diff(y(x), x))+(x^2-2)*y(x) = 0

and evaluate it where y(1) = 1 and y(2) = 2 to find y(3).


When I do the equation,


ODE := x^2*(diff(y(x), x, x))+x^3*(diff(y(x), x))+(x^2-2)*y(x) = 0

dsolve(ODE) y(x) = _C1/x+_C2*(-sqrt(Pi)*sqrt(2)*erf((1/2)*sqrt(2)*x)+2*x*exp(-(1/2)*x^2))/x


I have gotten an erf. Is this correct?

5 6 7 8 9 10 11 Last Page 7 of 20