Items tagged with ode ode Tagged Items Feed

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.

 

 

 

 

 

 

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)

 

Your remarks.

Thanks

 

 

 

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?

 

Hello,

I would like to solve the differential equation in the following link:

http://www.utdallas.edu/~frensley/technical/nanomes91/node2.html

 

Without using any explicit discretization. Is it possible to solve this equation with a Maple dsolve comand and some boundary condition option?

In ode solve command i generated a large array data. The output shows a large order matrix of this form

 

[110001x6 Matrix

Datatype:Anything

Storage:rectangular

order:Fortran_order].

 

I want to export this matrix into a notepad. Which can then be used for plotting in TecPlot. 

 

Looking for good response

 

 

Hi,

I have an ode like this

ODE:=(diff(T(x), x, x))+P*(S+a*(1-exp(-L*x))/L)*(diff(T(x), x))=0;

bcs:=T(0)=1,T(infinit)=0;

where P, S, L a, are all constants.

let assume that 

z=P/L^2*exp(-L*x);

subing z into the ode, we can have

ode1:=diff(T(z), z$2)+(1+z*a-P)*(diff(T(z), z)) = 0;

bcs1:=T(P/L^2)=1,T(0)=0;

Is it possible to find a closed form solution (T(x) in a compact form)?

hello,

restart:

ODE:=diff(T(z),z$2)+A1*(S-1/L+1/L*exp(-L*z))*diff(T(z),z)+A2*T(z)=0;

bcs:=T(0)=1,T(infinity)=0;

bcs:=T(0)=1,T(A3)=0;

dsolve({ODE,bcs});

where, A1, A2, A3, L, S are all constants.

i get an exact solution but is there any way around to get a more compact solution?

 

 

Hi,

I have a system of differential equations with boundary conditions:

diff(S(t), t) = -K(t)*S(t)/N, diff(K(t), t) = K(t)*S(t)/N, where S(T)=10, K(T)=20; If I would like to solve this system backward in time, is it right to re-write the system of original diff. equations in the following way:

diff(S(t), t) = K(t)*S(t)/N, diff(K(t), t) = -K(t)*S(t)/N and S(0)=10, K(0)=20 ( I simply changed the sign of the right hand of the equations and changed the boundary conditions to the inotoal ones).

Thanks,

Dmitry

 

 

Dear Experts,

When I run this code in maple I am facing with "Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging".

restart:
 
 unprotect('gamma');
 lambda:=5*10^5:
 mu:=0.003:
 beta:=4*10^(-10):
 delta:=0.2:
 alpha:=0.043:
 sigma:=alpha+delta:
 k:=6.24:
 gamma:=0.65:
 A[1]:=1:
 A[2]:=1:

ics := x[1](0)=1.7*10^8, x[2](0)=0,x[3](0)=400,psi[1](50)=0,psi[2](50)=0,psi[3](50)=0:

ode1:=diff(x[1](t), t)=lambda-mu*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)+delta*x[2](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)-sigma*x[2](t),
 diff(x[3](t), t) =(1+psi[3](t)*k*x[2](t)/A[2])*k*x[2](t)-gamma*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)*(-mu+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*k^2*x[2](t)-psi[1](t)*delta+psi[2](t)*sigma-psi[3](t)*(psi[3](t)*k^2/A[2]*x[2](t)+(1+psi[3](t)*k*x[2](t)/A[2])*k),
> 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)*gamma;

sol:=dsolve([ode1,ics],numeric, method = bvp[midrich]);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

Please help me to solve this equation on Maple.




Following my previous question

http://www.mapleprimes.com/questions/200627-Lssolve-Midpoint

I wrote the following code

 

restart:
Phiavg:=0.06;
lambda:=0.05;
Ha:=0;
NBT:=0.5;
Nr:=500;
#N[bt]:=cc*NBT+(1-cc)*4; ## cc between 0 and 1
N[bt]:=cc*NBT+(1-cc^2)*0.75;


                              0.06
                              0.05
                               0
                              0.5
                              500
                                           2
                    0.5 cc + 0.75 - 0.75 cc
eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(sigma-Nr*(phi(eta)-phi1[w])-(1-phi(eta))*T(eta)-Ha^2*u(eta))+((1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
eq2:=diff(T(eta),eta)-1/(k(eta)/k1[w]);
eq3:=diff(phi(eta),eta)-phi(eta)/(N[bt]*(1-gama1*T(eta))^2)*diff(T(eta),eta);
 /  d   /  d         \\      1                                 
 |----- |----- u(eta)|| + ------- (mu1[w] (sigma - 500 phi(eta)
 \ deta \ deta       //   mu(eta)                              

    + 500 phi1[w] - (1 - phi(eta)) T(eta)))

             /  d           \ /  d         \
      mu_phi |----- phi(eta)| |----- u(eta)|
             \ deta         / \ deta       /
    + --------------------------------------
                     mu(eta)                
                    /  d         \   k1[w]
                    |----- T(eta)| - ------
                    \ deta       /   k(eta)
                                       /  d         \            
                              phi(eta) |----- T(eta)|            
/  d           \                       \ deta       /            
|----- phi(eta)| - ----------------------------------------------
\ deta         /   /                       2\                   2
                   \0.5 cc + 0.75 - 0.75 cc / (1 - gama1 T(eta))
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):
rhop:=3880:
rhobf:=998.2:
cp:=773:
cbf:=4182:
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):
gama1:=0.00:
a[mu1]:=39.11:
b[mu1]:=533.9:
mu1[bf]:=9.93/10000:
a[k1]:=7.47:
b[k1]:=0:
k1[bf]:=0.597:
zet:=1:
phi1[w]:=phi0:
mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1);
eq2:=subs(phi(0)=phi0,eq2);
eq3:=subs(phi(0)=phi0,eq3);
/  d   /  d         \\   //                                    
|----- |----- u(eta)|| + \\0.0009930000000 + 0.03883623000 phi0
\ deta \ deta       //                                         

                      2\                                 
   + 0.5301627000 phi0 / (sigma - 500 phi(eta) + 500 phi0

                           \//               
   - (1 - phi(eta)) T(eta))/ \0.0009930000000

                                                   2\   
   + 0.03883623000 phi(eta) + 0.5301627000 phi(eta) / +

  /                                       /  d           \ /  d  
  |(0.03883623000 + 1.060325400 phi(eta)) |----- phi(eta)| |-----
  \                                       \ deta         / \ deta

         \\//                                        
   u(eta)|| \0.0009930000000 + 0.03883623000 phi(eta)
         //                                          

                          2\
   + 0.5301627000 phi(eta) /
           /  d         \     0.597 + 4.45959 phi0  
           |----- T(eta)| - ------------------------
           \ deta       /   0.597 + 4.45959 phi(eta)
                                        /  d         \
                            1. phi(eta) |----- T(eta)|
         /  d           \               \ deta       /
         |----- phi(eta)| - --------------------------
         \ deta         /                           2
                             0.5 cc + 0.75 - 0.75 cc  
Q:=proc(pp2,fi0) option remember; local res,F0,F1,F2,a,INT0,INT10,B;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if;
res := dsolve(subs(sigma=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}), numeric,output=listprocedure,initmesh=10, continuation=cc);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
INT0:=evalf(Int((abs(F0(eta)),eta=0..1)));
INT10:=evalf(Int(abs(F0(eta))*F1(eta),eta=0..1));
a[1]:=evalf(Int(F0(eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf),eta=0..1));
#a[1]:=evalf(Int((F0(eta),eta=0..1)));
a[2]:=(INT10/INT0-Phiavg)/Phiavg; #relative
[a[1],a[2]]
end proc:
Q1:=proc(pp2,fi0) Q(_passed)[1] end proc;
Q2:=proc(pp2,fi0) Q(_passed)[2] end proc;
proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
#Q(116,0.0041);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[130,0.01]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[43.55,0.39]);
tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5.65,0.00036]);
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[12,0.003]); # khoob ba 1
#tempe:=Optimization:-LSSolve([Q1,Q2],initialpoint=[5,0.01]);
                  HFloat(5.65), HFloat(3.6e-4)
           HFloat(5.650000070103341), HFloat(3.6e-4)
           HFloat(5.65), HFloat(3.600105456508193e-4)
     HFloat(29.63242379055208), HFloat(0.0205927592420527)
    HFloat(12.803902258015825), HFloat(0.006395385884750864)
    HFloat(12.803902403534572), HFloat(0.006395385884750864)
    HFloat(12.803902258015825), HFloat(0.00639539649402585)
   HFloat(12.804004931505949), HFloat(0.0063954867657199386)
    HFloat(12.804107604996073), HFloat(0.006395587646689013)
    HFloat(12.80400483062498), HFloat(0.006498160255844027)
    HFloat(12.803902157134855), HFloat(0.006498059374874952)
   HFloat(-1.0206939292143726), HFloat(-3.32764179807047e-4)
   HFloat(-1.0206939079125088), HFloat(-3.32764179807047e-4)
   HFloat(-1.0206939292143726), HFloat(-3.327536344433438e-4)
    HFloat(18.749500943683863), HFloat(0.01993840615828979)
    HFloat(3.9953780262640484), HFloat(0.00481041471606933)
     HFloat(6.166152606930136), HFloat(0.00703619658484674)
    HFloat(7.3193201827812295), HFloat(0.008218585352824569)
Error, (in Optimization:-LSSolve) complex value encountered
sigma:=tempe[2](1);
                          tempe[2](1)
phi0:=tempe[2](2);
                          tempe[2](2)
with(plots):

res2 := dsolve({eq1=0,eq2=0,eq3=0,u(1)=-lambda*D(u)(1),u(0)=lambda*D(u)(0),phi(0)=phi0,T(0)=0}, numeric,output=listprocedure,continuation=cc);
Error, (in dsolve/numeric/process_input) boundary conditions specified at too many points: {0, 1, 2}, can only solve two-point boundary value problems
G0,G1,G2:=op(subs(res2,[u(eta),phi(eta),T(eta)])):
ruu:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet)))/(Phiavg*rhop+(1-Phiavg)*rhobf);
phb:=evalf((Int(abs(G0(eta))*G1(eta),eta=0..1))) / evalf((Int(abs(G0(eta)),eta=0..1))) ;
TTb:=evalf(Int(abs(G0(eta))*G2(eta)*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1))/evalf(Int(abs(G0(eta))*(G1(eta)*rhop*cp+(1-G1(eta))*rhobf*cbf ),eta=0..1));
Error, invalid input: subs received res2, which is not valid for its 1st argument
                /  /1.                                        \
                | |                                           |
0.0008538922115 | |    |G0(eta)| (2881.8 G1(eta) + 998.2) deta|
                | |                                           |
                \/0.                                          /
                    /1.                       
                   |                          
                   |    |G0(eta)| G1(eta) deta
                   |                          
                  /0.                         
                  ----------------------------
                        /1.                   
                       |                      
                       |                      
                       |    |G0(eta)| deta    
                      /                       
                       0.                     
                                                              /Int(
                              1                               |     
------------------------------------------------------------- |     
  /1.                                                         |     
 |                                                            \     
 |              /             6                       6\            
 |    |G0(eta)| \-1.1752324 10  G1(eta) + 4.1744724 10 / deta       
/                                                                   
 0.                                                                 

                    /             6                       6\ , eta = 0. .. 1.)
  |G0(eta)| G2(eta) \-1.1752324 10  G1(eta) + 4.1744724 10 /                  

  \
  |
  |
  |
  /
#rhouu:=evalf((Int((G1(eta)*rhop+(1-G1(eta))*rhobf)*G0(eta),eta=0..1)));

odeplot(res2,[[eta,u(eta)/ruu],[eta,phi(eta)/phb],[eta,T(eta)/TTb]],0..1);
#odeplot(res2,[[eta,u(eta)],[eta,phi(eta)],[eta,T(eta)]],0..1);
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
rhou:=evalf((Int(abs(G0(eta))*(G1(eta)*rhop+(1-G1(eta))*rhobf ),eta=0..zet))):

Nub:=(1/G2(1))*(((1+a[k1]*abs(G1(0))+b[k1]*abs(G1(0))^2)/(1+a[k1]*Phiavg+b[k1]*Phiavg^2)));
                0.6905123602 (1 + 7.47 |G1(0)|)
                -------------------------------
                             G2(1)             
(rhs(res2(0.0000000000001)[3])-rhs(res2(0)[3]))/0.0000000000001;
Error, invalid input: rhs received res2(0.1e-12)[3], which is not valid for its 1st argument, expr
sigma;
                          tempe[2](1)
NBT;
                              0.5
>

 

the above code has been worked for NBT=0.6 and higher, whereas as NBT decreases, the code doesnt converge easily.

How can I fix this problem?

Thanks for your attention in advance

Amir

Dear Sirs,

I actually rigoruos to know what is the algorithm of BVP[midrich]? how it can obtain the solution of ODE with singularities?

 

Did anyone introduce a reference about the algorithm like this?

Thanks for your attention in advance

Amir

Is it possible that this expression has an elementary one (specifically the dilog's):

Y0:=(1/16)*(s*t*(exp(2*t)*s+exp(4*t)+1)*ln((exp(2*t)*s-(-s^2+1)^(1/2)+1)^16*(1+(-s^2+1)^(1/2))^16/((exp(2*t)*s+(-s^2+1)^(1/2)+1)^16*(1-(-s^2+1)^(1/2))^16))+s^3*t*(exp(4*t)+1)*ln((exp(2*t)*s-(-s^2+1)^(1/2)+1)^8*(1+(-s^2+1)^(1/2))^8/((exp(2*t)*s+(-s^2+1)^(1/2)+1)^8*(1-(-s^2+1)^(1/2))^8))+exp(2*t)*t*ln((exp(2*t)*s-(-s^2+1)^(1/2)+1)^32*(1+(-s^2+1)^(1/2))^32/((exp(2*t)*s+(-s^2+1)^(1/2)+1)^32*(1-(-s^2+1)^(1/2))^32))+4*((exp(4*t)+1)*s+2*exp(2*t))*(s^2+2)*dilog((-exp(2*t)*s+(-s^2+1)^(1/2)-1)/(-1+(-s^2+1)^(1/2)))-4*((exp(4*t)+1)*s+2*exp(2*t))*(s^2+2)*dilog((exp(2*t)*s+(-s^2+1)^(1/2)+1)/(1+(-s^2+1)^(1/2)))+((32*s^2*t+64*t)*exp(2*t)+16*(((t+1/8)*s^2+2*t+2)*exp(4*t)-(5/4)*s*exp(-2*t)-(1/8)*exp(-4*t)*s^2+(5/4)*s*exp(6*t)+(1/8)*s^2*exp(8*t)+(t-1/8)*s^2-2+2*t)*s)*arctanh((exp(2*t)-1)*(-1+s)/((-s^2+1)^(1/2)*(exp(2*t)+1)))+8*(-s^2+1)^(1/2)*((1/8)*s*(exp(4*t)+1)*ln((exp(4*t)*s+2*exp(2*t)+s)^12/s^12)+(1/8)*exp(2*t)*ln((exp(4*t)*s+2*exp(2*t)+s)^24/s^24)+(s^2-6*t-3)*exp(2*t)+((-(1/8)*s^2-3*t)*exp(4*t)+s*exp(-2*t)+(1/8)*exp(-4*t)*s^2+s*exp(6*t)+(1/8)*s^2*exp(8*t)-(1/8)*s^2-3*t)*s))/((s*exp(-2*t)+exp(2*t)*s+2)*(exp(4*t)*s+2*exp(2*t)+s)*((-s^2+1)^(1/2)+2*arctanh((-1+s)/(-s^2+1)^(1/2))))

Also I'm wondering since Y0 should solve the ode

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

with some constant C but I only get rubbish.

I ask this because I found that in another context this seems to be correct:

f1:=-(1/12)*Pi^2*((-s^2+1)^(1/2)-arccosh(1/s))/(-s^2+1)^(3/2)+(1/12)*arccosh(1/s)^3/(-s^2+1)^(3/2)-(1/4)*arccosh(1/s)^2/(-s^2+1)

f2:=(1/2)*((-s^2+1)^(1/2)*(polylog(2, s/(-1+(-s^2+1)^(1/2)))+polylog(2, -s/(1+(-s^2+1)^(1/2))))-polylog(3, s/(-1+(-s^2+1)^(1/2)))+polylog(3, -s/(1+(-s^2+1)^(1/2))))/(-s^2+1)^(3/2)

and f1=f2

but maple doesnt convert it.

Also maple has trouble to convert

2*arctanh(sqrt((1-s)/(1+s)))=arccosh(1/s)

everywhere: 0<s<1

Hi,

I get the error in the following code

restart:

gama1:=0.01:

zet:=0;
#phi0:=0.00789:
Phiavg:=0.02;
lambda:=0.01;
Ha:=1;


                               0
                              0.02
                              0.01
                               1
rhocu:=2/(1-zet^2)*int((1-eta)*rho(eta)*c(eta)*u(eta),eta=0..1-zet):

eq1:=diff(u(eta),eta,eta)+1/(mu(eta)/mu1[w])*(1-Ha^2*u(eta))+((1/(eta)+1/mu(eta)*(mu_phi*diff(phi(eta),eta)))*diff(u(eta),eta));
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]/(eta)*diff(T(eta),eta) ));
eq3:=diff(phi(eta),eta)+phi(eta)/(N[bt]*(1+gama1*T(eta))^2)*diff(T(eta),eta);
      /  d   /  d         \\   mu1[w] (1 - u(eta))
      |----- |----- u(eta)|| + -------------------
      \ deta \ deta       //         mu(eta)      

           /             /  d           \\               
           |      mu_phi |----- phi(eta)||               
           | 1           \ deta         /| /  d         \
         + |--- + -----------------------| |----- u(eta)|
           \eta           mu(eta)        / \ deta       /
                                /      /                        
                                |      |                        
/  d   /  d         \\     1    |      |  rho(eta) c(eta) u(eta)
|----- |----- T(eta)|| + ------ |k1[w] |- ----------------------
\ deta \ deta       //   k(eta) |      |         5000 p2        
                                \      \                        

                                /  d           \
     (a[k1] + 2 b[k1] phi(eta)) |----- phi(eta)|
                                \ deta         /
   + -------------------------------------------
                                          2     
         1 + a[k1] phi1[w] + b[k1] phi1[w]      

            /  d         \\\
     k(eta) |----- T(eta)|||
            \ deta       /||
   + ---------------------||
           k1[w] eta      ||
                          //
                                      /  d         \
                             phi(eta) |----- T(eta)|
          /  d           \            \ deta       /
          |----- phi(eta)| + ------------------------
          \ deta         /                          2
                             N[bt] (1 + 0.01 T(eta))
mu:=unapply(mu1[bf]*(1+a[mu1]*phi(eta)+b[mu1]*phi(eta)^2),eta):
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta):
rhop:=3880:
rhobf:=998.2:
cp:=773:
cbf:=4182:
rho:=unapply(  phi(eta)*rhop+(1-phi(eta))*rhobf ,eta):
c:=unapply(  (phi(eta)*rhop*cp+(1-phi(eta))*rhobf*cbf )/rho(eta) ,eta):
mu_phi:=mu1[bf]*(a[mu1]+2*b[mu1]*phi(eta)):

a[mu1]:=39.11:
b[mu1]:=533.9:
mu1[bf]:=9.93/10000:
a[k1]:=7.47:
b[k1]:=0:
k1[bf]:=0.597:
zet:=0.5:
#phi(0):=1:
#u(0):=0:
phi1[w]:=phi0:
N[bt]:=0.2:
mu1[w]:=mu(0):
k1[w]:=k(0):

eq1:=subs(phi(0)=phi0,eq1):
eq2:=subs(phi(0)=phi0,eq2):
eq3:=subs(phi(0)=phi0,eq3):

#A somewhat speedier version uses the fact that you really need only compute 2 integrals not 3, since one of the integrals can be written as a linear combination of the other 2:
Q:=proc(pp2,fi0) local res,F0,F1,F2,a,INT0,INT10,B;
global Q1,Q2;
print(pp2,fi0);
if not type([pp2,fi0],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve(subs(p2=pp2,phi0=fi0,{eq1=0,eq2=0,eq3=0,u(1)=lambda/(phi(1)*rhop/rhobf+(1-phi(1)))*D(u)(1),D(u)(0)=0,phi(1)=phi0,T(1)=0,D(T)(1)=1}), numeric,output=listprocedure):
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)])):
INT0:=evalf(Int((1-eta)*F0(eta),eta=0..1-zet));
INT10:=evalf(Int((1-eta)*F0(eta)*F1(eta),eta=0..1-zet));
B:=(-cbf*rhobf+cp*rhop)*INT10+ rhobf*cbf*INT0;
a[1]:=2/(1-zet^2)*B-10000*pp2;
a[2]:=INT10/INT0-Phiavg;
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;
#The result agrees very well with the fsolve result.
#Now I did use a better initial point. But if I start with the same as in fsolve I get the same result in just about 2 minutes, i.e. more than 20 times as fast as fsolve:

Q1:=proc(pp2,fi0) Q[1](_passed) end proc;
Q2:=proc(pp2,fi0) Q[2](_passed) end proc;
Optimization:-LSSolve([Q1,Q2],initialpoint=[6.5,exp(-1/N[bt])]);


proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
proc(pp2, fi0)  ...  end;
              HFloat(6.5), HFloat(0.006737946999)

 

 

the error is :

Error, (in Optimization:-LSSolve) system is singular at left endpoint, use midpoint method instead

how can I fix it.

Thanks

 

Amir

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