# Items tagged with odeode Tagged Items Feed

### 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 ;

### ode plot help??...

March 12 2014
0 5

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?

### convert 2nd order ODE to system of First ODE ...

March 09 2014
1 2

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

### Numerica Dsolve for Boundary values problem...

March 09 2014
1 8

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 ### Billiards with Maple... February 24 2014 0 15 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. ### Adams–Moulton methods to solve ODE... February 23 2014 0 3 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. ### Midpoint Method using for solving Ordinary differe... February 22 2014 1 16 Dear all; Thanks ifor looking and help me in my work. Your remarks are welcome.Description: This routine uses the midpoint method to approximate the solution of the differential equation$y'=f(x,y)$with the initial condition$y = y[0]$at$x = a$and recursion starting value$y = y[1]$at$x = a+h$. The values are returned in$y[n]$, the value of$y$evaluated at$x = a + nh$. Arguments: \begin{itemize} \item$f$the integrand, a function of a two variables \item$y[]$On input$y[0]$is the initial value of$y$at$x = a$, and$y[1]$is the value of$y$at$x = a + h$, \item on output for$i \geqslant 2\$
$$y[i] = y[i-2] + 2h f(x[i],y[i]); \quad \quad x[i] = a + i h.$$
\end{itemize}

CODE USING MAPLE

Midpoint-Method=proc(f,a,b, N)

h:=(b-a)/N;
x[0]:=a;
y[0]:=1:

for n from 2 to N do
x[n] := a+n*h;
y[n+1] = y[n-1] +  2h f( x[n], y[n] );
od:
// Generate the sequence of approximations for the Improved Euler method
data_midpoint := [seq([x[n],y[n]],n=0..N)]:
//write the function;
F:=(t,y)-> value of function ;

//Generate plot which is not displayed but instead stored under the name out_fig for example
out_fig := plot(data_midpoint,style=point,color=blue)

Thanks

### Using numerical sol of ODE in another ODE...

February 19 2014
0 2

Say we solve numerically and ODE using Maple. Say ode1:= { diff(Q(x),x)= Q(x)/3x , Q(1)=1 }

The solution is a procedure so now suppose we have another ODE where the solution appears.Say  ode2:= { diff(f(x),x)= Q(x)*x , f(1)=4}.

To extract the solution of the first ODE I set sol1:=dsolve(ode1,numeric) and Q:=proc(x) local s: return rhs(sol(x)[2]): end proc:

But now I got an error message when I trying sol2:=dsolve(ode2,numeric).

Is it possible to use a procedure in the definition of the ODE one wants to solve?

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