Items tagged with ode ode Tagged Items Feed

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.  

 Math_4710_Hilbert_16.mw

Hi all

I dont know why some Z1 appears on the screen and the code does not converge.
please help me
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"]);
plot(pade(diff(f(tau),tau), tau, [7, 7]),tau=0..5,color=blue,style=point,symbol=circle,symbolsize=7,labels=["tau","velocity"]);

 

 

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.

restart:

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

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

dsolve({ODE,bcs});

             u(y) = 0

 

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

 

 

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;
RKadaptivestepsize((x,y)->2*x-y,0,1, epsilon,20)

Hi.

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

Thank you.

Dear all,

Thank you for your Help.

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.

  epsilon:=(x0,h)->sqrt( 1/(N+1)*add( (Heun1(f,x0,i)-Heun2(f,x0,i))^2,i=0..N ) );

  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;

 

 

Thanks you for your help.


                                

                        

 

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.

Thanks for your attentions in advance

Amir

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 ;

 

 

ode plot help??...

March 12 2014 J4James 175

restart:

Eq1:=diff(f(x),x)=1;

f:=1:

plot(1,x=0..10);

cs:=(f)(0)=0;

p:=dsolve({Eq1,cs},numeric);

odeplot(p, [x,diff(f(x),x)], 0..10);

why I cannot get a plot diff(f(x),x) like this?

 

whats going on here?

 

 

 

Dear all;

Please Have some one an idea to transform or convert 2nd order ODE to system of First ODE ( of course using maple).

Thanks

 

Hi,

I wrote the following code which is properly run

 


restart:

# parametrs

MUR:=(1-phi)^2.5:
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:

dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:

k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:

phi:=0.00:
binfinitive:=6: Pr:=7: lambda:=0:


with(plots):
pppe:=dsolve( {dqu1=0,dqu2=0,dqu3=0,T(0)=1,T(binfinitive)=0,f(0)=0,D(f)(0)=lambda,D(f)(binfinitive)=0,h(binfinitive)=0}, numeric );
-pppe(0);
print(odeplot(pppe,[x,diff(f(x),x)],0..binfinitive,color=black,numpoints=400));
print(odeplot(pppe,[[x,diff(f(x),x)]],0..binfinitive,color=black,numpoints=400));
print(odeplot(pppe,[[x,T(x)]],0..binfinitive,color=black,numpoints=400));


However, in some range of parameters, I must increase the value of binfinitive (for example binfinitive=50). however, my code is doesnt converge for higher values of 10 (at most). Can anyone change this algorithm in a way that it insensitive to the value of binfinitive?

Many thanks for your attention in advance

 

Amir

Is it possible to describe billiards in the unit square (see http://en.wikipedia.org/wiki/Dynamical_billiards for info) in terms of events and if? For example, let the dynamical system
{x''(t)=0,y''(t)=0, x(0)=0.5, x'(0)=0.2, y(0)=0.5, y'(0)= sqrt(3)/5}
be given, where the trajectory [x(t), y(t)] satisfies the usual laws of reflection http://en.wikipedia.org/wiki/Reflection_%28physics%29 when reaching the boundary {x=0}, {x=1}, {y=0}, {y=1}. The same question about billiards in the unit disk.

Dear all

Is there any one can help me to find  the Maple code to solve ODE : y'(x)=f(x,y(x))  using n-step  Adams-Moulton Methods.

The code exist  with mathematica in this link:

http://mathfaculty.fullerton.edu/mathews/n2003/AdamsBashforthMod.html

there is also the code of this method with Matlab, see please:

 http://www.math.mcgill.ca/gantumur/math579w10/matlab/abm4.m

 

But I want file.mw ( with maple)

Thank you very much for helping me.

 

 

 

 

 

 

1 2 3 4 5 6 7 Last Page 3 of 18