# Items tagged with odeode Tagged Items Feed

### Finding zeroes of a numerical solution of an ODE i...

April 22 2014
2 4

I have numerically solved a system of ODEs and plotted the graphs of a[j](t) for each j=0..21.

It was clear from the picture that each a[j] has a unique zero. Is there a maple command to

locate these zeroes?

### ODE with constraint...

April 20 2014
0 1

Hi

I have the following ODEs

Restart:
a:=0.13:
b:=0.41:
reynolds:=1.125*10^8:
Eq1:=diff(f(x),x\$3)+diff(f(x),x\$2)*f(x)+b^2*sqrt(2*reynolds)*diff(diff(f(x),x\$2)*f(x)*x^2,x\$1);
Eq2:=diff(g(x),x\$3)+diff(g(x),x\$2)*g(x)+c*a^2*sqrt(2*reynolds)*diff(diff(g(x),x\$2)*x,x\$1);
f(0)=0;
D(f(0)):=0;
# continuity condition
g(0.1*c)=f(0.1*c):
D(g(0.1*c))=D(f(0.1*c)):
(D@2)(g(0.1*c))=(D@2)(f(0.1*c)):

the value of c is unknown which must be obtained via D(g(c)):=1;

How can I solve it?

Amir

### Error, (in dsolve/numeric/bvp/convertsys) too few ...

April 19 2014
0 6

im solving 4 ODe with boundary conditions.. i got this error Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 8, got 7

### How to plot a graph...

April 16 2014
1 4

i solve 4 ODE with boundary condition.. i try to plot a graph F(eta) with different value of M.. but it doesnt comeout.. anyone can help me please??

 (1)

 (2)

 (3)

 (4)

 (5)

 (6)

 (7)

 (8)

 (9)

 (10)

 (11)

### Error, (in dsolve/numeric/bvp/convertsys) unable t...

April 08 2014
0 13

i am trying to solve 6 ODE with boundary condition

 (1)

 (2)

 (3)

 (4)

 (5)

 (6)

 (7)

 (8)

 (9)

 (10)

 (11)

 (12)

 (13)

 (14)

 (15)

 (16)

 (17)

 (18)

 (19)

 (20)

 (1)

 (2)

 (3)

 (4)

 (5)

 (6)

 (7)

 (8)

 (9)

 (10)

 (11)

 (12)

 (13)

 (14)

 (15)

 (16)

 (17)

 (18)

 (19)

 (20)

then i get this error

Error, (in dsolve/numeric/bvp/convertsys) unable to convert to an explicit first-order system

i dont know where i need to change after view it one by one..

### Problems with limit cycles and maple 16 ...

April 04 2014
0 1

I am having trouble printing out a limit cylce on maple 16.  I have the attached file and if anybody could look at it and perhaps help me out it would be greatly appreciated.  The first limit cycle is supposed to look somewhat like the second one.  I'v tried many different things but nothing seems to be working.  an explenation would also be nice too.  if the file does not open correctly also let me know. thank you very much.

### problem with convergence ...

March 30 2014
0 13

Hi all

I dont know why some Z1 appears on the screen and the code does not converge.
thanks alooooot

restart;
n:=3;
nn:=3;
m:=1;
BB:=1;
BINF:=5:
pr:=7;
digits:=10;
>
eq1:=diff(f(tau),tau\$3)+((3/5)*f(tau)*diff(f(tau),tau\$2))-(1/5)*(diff(f(tau),tau\$1))^2+((2/5)*tau*diff(h(tau),tau\$1))-((2/5)*h(tau))-BB*diff(f(tau),tau\$1)=0;
eq11:=(1/pr)*diff(h(tau),tau\$3)+(3/5)*f(tau)*diff(h(tau),tau\$2)=0;

h(tau):=sum(p^i*h[i](tau),i=0..nn);
f(tau):=sum(p^i*f[i](tau),i=0..n);

H1:= p*(diff(f(tau),tau\$3)+((3/5)*f(tau)*diff(f(tau),tau\$2))-(1/5)*(diff(f(tau),tau\$1))^2+((2/5)*tau*diff(h(tau),tau\$1))-((2/5)*h(tau))-BB*diff(f(tau),tau\$1))+(1-p)*(diff(f(tau),tau\$3)):
H11:= p*((1/pr)*diff(h(tau),tau\$3)+(3/5)*f(tau)*diff(h(tau),tau\$2))+(1-p)*(diff(h(tau),tau\$3)):
>
eq2:=simplify(H1):
eq22:=simplify(H11):
eq3:=collect(expand(eq2),p):
eq33:=collect(expand(eq22),p):
eq4:=
convert(series(collect(expand(eq2), p), p, n+1), 'polynom');
eq44:=
convert(series(collect(expand(eq22), p), p, n+1), 'polynom');
for i to n do
s[i] := coeff(eq4, p^i) ;
print (i);
end do:
for i to nn do
ss[i] := coeff(eq44, p^i) ;
print (i);
end do:
s[0]:=diff(f[0](tau), tau\$3);
ss[0]:=diff(h[0](tau), tau\$3);
ics[0]:=f[0](0)=0, D(f[0])(0)=0, D(f[0])(BINF)=0;
icss[0]:=h[0](BINF)=0, D(h[0])(0)=1, D(h[0])(BINF)=0;

dsolve({s[0], ics[0]});
f[0](tau):= rhs(%);
#dsolve({ss[0], icss[0]});
h[0](tau):= -exp(-tau); #;rhs(%);

>
>
for i to n do
f[ii-1](tau):=convert(series(f[ii-1](tau), tau, nn+1), 'polynom');
h[ii-1](tau):=convert(series(h[ii-1](tau), tau, nn+1), 'polynom');
s[i]:=simplify((s[i]));
ics[i]:=f[i](0)=0, D(f[i])(0)=0, D(f[i])(BINF)=0;
dsolve({s[i], ics[i]});
f[i](tau):=rhs(%);
ss[i]:=(ss[i]);
icss[i]:=h[i](BINF)=0, D(h[i])(0)=0, D(h[i])(BINF)=0;
dsolve({ss[i], icss[i]});
h[i](tau):=rhs(%);

end do;

f(tau):=sum((f[j])(tau),j=0..n);
with(numapprox):

plot(diff(f(tau),tau),tau=0..5,color=blue,style=point,symbol=circle,symbolsize=7,labels=["tau","velocity"]);

### Second order direction field problems...

March 29 2014
0 1

I'm trying to plot the direction field of the second order differential equation x''=x'-cos(x) using dfieldplot:

> with(DEtools); with(plots);
> f1 := (x, y) options operator, arrow; diff(x(t), t)-cos(x(t)) end proc;
/ d \
(x, y) -> |--- x(t)| - cos(x(t))
\ dt /
> dfieldplot([diff(x(t), t) = y(t), diff(y(t), t) = f1(x(t), y(t))], [x(t), y(t)], t = -2 .. 2, x = -2 .. 2, y = -2 .. 2);
Error, (in DEtools/dfieldplot) cannot produce plot, non-autonomous DE(s) require initial conditions.
>

The error I'm getting says I need initial conditions, but I wasn't provided with any. Is there another way to plot this? Sorry if this is dumb question, but I've only ever plotted first order equations.

### Does this ode has a non-trivial sol?...

March 29 2014
1 3

restart:

ODE:=diff((-diff(u(y),y))^n,y)=A;

bcs:=D(u)(0)=0,u(h)=0;

dsolve({ODE,bcs});

u(y) = 0

### Solving ODE's system and couldnt get any result. ...

March 23 2014
2 21

Hello everyone,

I have a problem solving with ODE's system solving. I have 2 equation and 4 initial  conditions. When i calculate like that,u can look this file.  beamsolvingfortim.mw. It is working it is giving me T1 and T2 equations depends on time.

In 2. system which i have a problem i want to calculate this equations depends on x(displacement). I have again 2 equation and 4 boundry conditions. it is solving the ODE'S system without boundry conditions. (It is giving with C1 C2 C3 and I)Problem is when i want to find its values with boundry conditions it is not giving a result. Is there a problem with complex number(I) or boundry conditions are not enough?

> evalf(dsolve({sys}));

{X1s(x) = -3.060206320 + _C1 exp(-0.3487988669 x)

+ _C2 exp(0.3487988669 x) + _C3 exp(-0.3563227426 I x), X2s(x) =

-0.6321326989 _C1 exp(-0.3487988669 x)

+ 0.6321326989 _C2 exp(0.3487988669 x)

- 0.6053448484 I _C3 exp(-0.3563227426 I x)}

When i write like that, nothing is happening. I also upload the files. If u can help , i would be really appreciate.

eq2 := dsolve({A, B, bc}, [X1s(x), X2s(x)])

beamsolving(thick.mw

### how can I use if when I solve this ode...

March 20 2014
1 12

Dear all;

Please I need your help to find the error in my code.

I want to solve an ode, with condition on step size.

ode := diff(y(x), x) = 2*x+y(x);
f:=(x,y)->2*x-y;

analyticsol := rhs(dsolve({ode, y(0) = 1}));
RKadaptivestepsize := proc (f, a, b, epsilon, N)
local x, y, n, h,k,z,R,p;
p:=2;
h := evalf(b-a)/N; ## we begin with this setpsize
x[0] := a; y[0] := 1; ## Initialisation
for n from 0 to N-1 do  ##loop
x[n+1] := a+(n+1)*h;  ## noeuds
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h, y[n]+h*k[1]);
k[3] := f(x[n]+h/2, y[n]+h/4*(k[1]+k[2]));
z[n+1] := z[n]+(h/2)*(k[1]+k[2]);## 2-stage runge Kutta.
y[n+1] := y[n]+(h/6)*(k[1]+k[2]+4*k[3]);
R:=abs(y[n+1]-z[n+1]); ## local erreur
hstar:=sqrt(epsilon/R)
if R=<=epsilon    then
x[n] := x[n+1]+h;
y[n]:=y[n+1];
n:=n+1;

else

h:=hstar;
end if
end do;
[seq([x[n], y[n]], n = 0 .. N)];
[seq([x[n], z[n]], n = 0 .. N)];
end proc:

epsilon:=1e-8;

### adaptative step size runge Kutta...

March 20 2014
0 3

Hi.

Please, I need a code in maple for adaptative setp size control for runge Kutta.

Thank you.

### Plot Error of ODE. ...

March 14 2014
0 5

Dear all,

h: stepsize;

x in [0,x0];

I give all the step of my code, but I think there is a mistake. I wait for your Help.

I would like to compute the error between  Method Huen with step size h and step size 2h using the definition of epsilon given below:

## The error written epsilon(x0,h)= sqrt(1/(N+1) * sum_i=0^N  (y_i^{2h}-y_(2i)^h)^2 ), where y_i^(2h) is the approximation of y(i*2*h).

## We want : loglog epsilon versus h.

f:=(x,y)=1/(1+cos(y));

ode:=diff(y(x),x)=f(x,y);

ic:=y(0)=1;  h:=x0/(2*N);

## Method Heun with step size 2h

> Heun1 := proc (f, x0,)

local x, y, i, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf((1/2)*x0/N);

x[0] := 0;

y[0] := 1;

for i from 0 to N do

x[i+1] := (2*i+2)*h;

k[1] := f(x[i], y[i]);

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+h*((1/2)*k[1]+(1/2)*k[2]);

end do;

[seq([x[i], y[i]], i = 0 .. N)];

end proc;

### Now Heun with step size h  ( the same h)

> Heun2 := proc (f, x0,)

local x, y, i, h, k;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf((1/2)*x0/N);

x[0] := 0;

y[0] := 1;

for i from 0 to N do

x[i+1] := (i+1)*h;

k[1] := f(x[i], y[i]);

k[2] := f(x[i]+h, y[i]+h*k[1]);

y[i+1] := y[i]+h*((1/2)*k[1]+(1/2)*k[2]);

end do;

[seq([x[2*i], y[2*i]], i = 0 .. N)];

end proc;

### widen the region of convergence ...

March 14 2014
0 2

Dear collegues

I wrote the following code

restart:
Digits := 15;
a[k]:=0;
b[k]:=7.47;
a[mu]:=39.11;
b[mu]:=533.9;
mu[bf]:=9.93/10000;
k[bf]:=0.597;
ro[p]:=3880 ;
ro[bf]:= 998.2;
c[p]:= 773;
c[bf]:= 4182;
#mu[bf]:=1;
Gr[phi]:=0; Gr[T]:=0;
#dp:=0.1;
Ree:=1;
Pr:=1;
Nbt:=cc*NBTT+(1-cc^2)*6;

#######################
slip:=0.1;         ####
NBTT:=2;           ####
lambda:=0.1;       ####
phi_avg:=0.02;    ####
#######################

eq1:=diff( (1+a[mu]*phi(eta)+b[mu]*phi(eta)^2)*diff(u(eta),eta),eta)+dp/mu[bf]+Gr[T]*T(eta)-Gr[phi]*phi(eta);
eq2:=diff((1+a[k]*phi(eta)+b[k]*phi(eta)^2)*diff(T(eta),eta),eta)+lambda*T(eta)/k[bf];
eq3:=diff(phi(eta),eta)+1/Nbt*diff(T(eta),eta);
Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10;
global Q1,Q2;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({subs(dp=pp2,eq1)=0,eq2=0,eq3=0,u(0)=slip*D(u)(0),u(1)=-slip*D(u)(1),D(T)(0)=0,D(T)(1)=1,phi(0)=fi0}, numeric,output=listprocedure,continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
INT0:=evalf(Int(F0(eta),eta=0..1));
INT10:=evalf(Int(F0(eta)*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta),eta=0..1))-Ree*Pr;;
a[2]:=INT10/INT0-phi_avg;
Q1(_passed):=a[1];
Q2(_passed):=a[2];
if type(procname,indexed) then a[op(procname)] else a[1],a[2] end if
end proc;
Q1:=proc(pp2,fi0) Q[1](_passed) end proc;
Q2:=proc(pp2,fi0) Q[2](_passed) end proc;
Optimization:-LSSolve([Q1,Q2],initialpoint=[0.3,0.0007]);

se:=%[2];
res2 := dsolve({subs(dp=se[1],eq1)=0,eq2=0,eq3=0,u(0)=slip*D(u)(0),u(1)=-slip*D(u)(1),D(T)(0)=0,D(T)(1)=1,phi(0)=se[2]}, numeric,output=listprocedure,continuation=cc);
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
TTb:=evalf(Int(G0(eta)*G2(eta)*(G1(eta)*ro[p]*c[p]+(1-G1(eta))*ro[bf]*c[bf] ),eta=0..1))/evalf(Int(G0(eta)*(G1(eta)*ro[p]*c[p]+(1-G1(eta))*ro[bf]*c[bf] ),eta=0..1));
with(plots):
odeplot(res2,[[eta,phi(eta)/phi_avg]],0..1);
odeplot(res2,[[eta,T(eta)/TTb]],0..1);
odeplot(res2,[[eta,u(eta)/(Ree*Pr)]],0..1);

res2(1);
Nuu:=(1/TTb);
1/((1+a[k]*G1(1)+b[k]*G1(1)^2)/(1+a[k]*phi_avg+b[k]*phi_avg^2));
(1/TTb)*(((1+a[k]*G1(1)+b[k]*G1(1)^2)/(1+a[k]*phi_avg+b[k]*phi_avg^2)));
>

I want to run the code for the value of NBTT in the range of 0.2 to 10. this code gave the results in the range of 4-10 easily. So, I used the continuation which improve the range of the results between 2-10. However, I coudnt gave the results when 0.2<NBTT<2. Would you please help me in this situation.

Also, It is to be said that the values of phi should be positive. in some ranges, I can see that phi(1) is negative. Can I place a condition in which the values phi restricted to be positive.

Amir

### Loop In MAPLE and test of stop...

March 13 2014
0 2

Dear all;

Than you for help.

how  many steps are required to achieve a error of 1.e-3 in the numerical value of y(1).

Here The 3 -step procedure  Range Kutta Method.

## Exact  solution

### We will modifty N ( number of steps to get error =10^(-3). )

## Procedure Range Kutta

> RK3 := proc (f, a, b, y0, N)

local x, y, n, h, k, vectRK3;

y := Array(0 .. N);

x := Array(0 .. N);

h := evalf(b-a)/N;

x[0] := a; y[0] := 1;

for n from 0 to N-1 do

x[n+1] := a+(n+1)*h;

k[1] := f(x[n], y[n]);

k[2] := f(x[n]+(1/2)*h, y[n]+(1/2)*h*k[1]);

k[3] := f(x[n]+h, y[n]+h*(-k[1]+2*k[2]));

y[n+1] := y[n]+(1/6)*h*(k[1]+4*k[2]+k[3])

end do;

[seq([x[n], y[n]], n = 0 .. N)]; y[1];

end proc;

## Now  we compute the error between y(1) and exact  solution for different value of  N

### I have a problem in this part

errorRk3 := array(1 .. 29);
for N from  2 to 30 do

errorrRk3[N] := abs(eval(rhs(res), x = 1)-RK3((x,y)->-y,0,1,N));

if errorrRk3[N] =10^{-3} end ;
end  do ;

 5 6 7 8 9 10 11 Last Page 7 of 22
﻿