Items tagged with numeric numeric Tagged Items Feed

Let us consider the definite integral

J:=int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1);

Maple fails with it, Mathematica 10.1 finds it in terms of  special functions. Let us look at the integrand:
plot(x-(-x^5+1)^(1/5), x = 0 .. 1);

We see the expression under the modulus changes its sign at the unique point of RealRange(0,1). Therefore

solve(x-(-x^5+1)^(1/5));


Then

J:= int(-x+(-x^5+1)^(1/5), x = 0 .. (1/2)*2^(4/5))+int(x-(-x^5+1)^(1/5), x = (1/2)*2^(4/5) .. 1);

which outputs a complicated expression

(1/8)*2^(4/5)*(4*hypergeom([-1/5, 1/5], [6/5], 1/2)-2^(4/5))+(1/2)*2^(4/5)*((1/2)*2^(1/5)-(1/4)*2^(4/5))-(1/25)*Pi*csc((1/5)*Pi)*(-(25/2)*sin((1/5)*Pi)*GAMMA(4/5)*2^(4/5)*hypergeom([-1/5, 1/5], [6/5], 1/2)/Pi+(5/4)*sec((3/10)*Pi)*cos((1/10)*Pi)*2^(3/5)*Pi^(1/2)*csc((3/10)*Pi)/GAMMA(7/10))/GAMMA(4/5).

At the same time we have

int(abs(x-(-x^5+1)^(1/5)), x = 0 .. 1, numeric);

                          0.5000000000

How to obtain 1/2 symbolically?






Hi guys,

I'm studying a system of six differential equations. Given the fact that the system cannot be solved symbolically, I've tried the numeric procedure, and it works. I proceeded like this :

soleqd:=dsolve(sysd2,numeric,var);

then i checked if maple could calculate the solutions for given values of t. It works for t=0, t=0.5, t=1,t=2,...,t=5. The solutions are all real numbers.

But when i try to draw a graphic representation of the solutions, it doesn't work. I do :

ff1:=t->subs(soleqd(t),u[1](t));
gg1:=t->subs(soleqd(t),nu[1](t));

Then :

plot(['ff1(t)','gg1(t)',t=0..5],u[1]=0...1,nu[1]=0...1);

(The square brackets are indices)

Now maple answers that it is "unable to evaluate the function to numeric values in the region". I went to the help page but no solution seems to work. I can't figure it out by myself. Does anybody notice something wrong with my code ?

Thank you for your time,

Best regards,

Louis

restart;

with(plots);

Eq1 := diff(f(eta), eta, eta, eta)+f(eta)*(diff(f(eta), eta, eta))-(diff(f(eta), eta))^2-M^2*(diff(f(eta), eta))+B(f(eta)*(diff(f(eta), eta, eta))*(diff(f(eta), eta))-f(eta)^2*(diff(f(eta), eta, eta, eta))) = 0;

Eq2 := (diff(theta(eta), eta, eta))/Pr+f(eta)*(diff(theta(eta), eta))-2*(diff(f(eta), eta))*theta(eta) = 0;

Pr := 1

M := 1

S := 0

epsilon := 1

blt := 10

bcs1 := f(0) = S, (D(f))(0) = epsilon, (D(f))(blt) = 0;

bcs2 := theta(0) = 1, theta(blt) = 0;

L := [0, .2, .4, .6, .8, 1.2];

for k to 6 do R := dsolve(eval({Eq1, Eq2, bcs1, bcs2}, B = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure); X1 || k := rhs(R[3]); X2 || k := rhs(R[4]); Y1 || k := rhs(R[5]); Y2 || k := -rhs(R[6]) end do:

print([(X2 || (1 .. 6))(0)])

Consider the equation

eqn1:=5^(2*abs(y-1)+2) = (1/15625)*625^abs(y-1);

and note that

eval( eqn1, y=-3);
eval( eqn1, y=5);

both return valid solutions. So why does

testeq(eqn1);

throw a numeric exception.

Given the manual description of testeq(), I would expect it to return either true or FAIL - but not a numeric exceprion!!!!

Any ideas?

 

Please, i need help USING ODE1 and ODE2 with given BCS and Pr=0.714

it is needed to to generates   

                   [-0.2], [0.51553], [0.4000]
                  [-0.1], [0.57000], [0.4371]
                   [0.], [0.62756], [0.4764]
                  [0.1], [0.68811], [0.5176]
                  [0.2], [0.75153], [0.5609]

but it is generarting

                  [-0.2], [0.51553], [0.42342]
                  [-0.1], [0.57000], [0.46114]
                   [0.], [0.62756], [0.50088]
                  [0.1], [0.68811], [0.54261]
                  [0.2], [0.75153], [0.58628]

the values of D(theata)(0) is wrong. Please i need HELP. this the code below that i use:
>restart;
>with (plots):ode1:=diff(f(eta),eta,eta,eta)+f(eta)*diff(f(eta),eta,eta)-M*diff(f(eta),eta)=0:

>ode2:=diff(theta(eta),eta,eta)+Pr*f(eta)*diff(theta(eta),eta)=0:

>bcs1:= f(0)=w,D(f)(0)=1,D(f)(10)=0:
>bcs2theta(10)=0,theta(0)=1:
>fixedparameter1:=[M=0.0]:
>ode3:=eval(ode1,fixedparameter1):
>fixedparameter2:=[Pr=0.714]:
>ode4:=eval(ode2,fixedparameter2):

>G:=[-0.2,-0.1,0.0,0.1,0.2]:
>for ode3 and ode4:
  for k from 1 to 5 do
  sol_All:=dsolve(eval({ode3,ode4,bcs1},w=G[k]),    [f(eta),theta(eta)],numeric,output=listprocedure);
Y_sol||k:= -rhs(sol_All[4]);
YP_sol||k:=-rhs(sol_All[6]);
end do:
>Digits:=5:

>for k from 1 to 5 do
evalf([G[k]]),evalf([(Y_sol||k(0))]),evalf([YP_sol||k(0)]);    
od;
                  [-0.2], [0.51553], [0.42342]
                  [-0.1], [0.57000], [0.46114]
                   [0.], [0.62756], [0.50088]
                  [0.1], [0.68811], [0.54261]
                  [0.2], [0.75153], [0.58628]

restart;
pp:=-55471918776960000*tanh((1/3220)*sqrt(10368400-cp^2)*Pi*x/cp)+5350094400*tanh((1/3220)*sqrt(10368400-cp^2)*Pi*x/cp)*cp^2-129*tanh((1/3220)*sqrt(10368400-cp^2)*Pi*x/cp)*cp^4+2670899840*tanh((1/6450)*sqrt(41602500-cp^2)*Pi*x/cp)*sqrt(41602500-cp^2)*sqrt(10368400-cp^2);
Student[Calculus1]:-Roots(subs(x=8000,pp),cp=1..3220,numeric);
p1:=proc(v)
option hfloat;
local a;
a:=Student[Calculus1]:-Roots(subs(x=v,pp),cp=1..3220,numeric);
if nops(a)>=1 then seq([v,a[i]],i=1..nops(a));
end if;
end proc:
SS1:=[seq(p1(i),i=3500..20000,200)]:
plot(SS1,style=point,gridlines);

The final figure is different between maple12 and maple17.

On 17, unwanted points apprear.

is it a bug?

Dear Colleges

I have a problem with the following code. As you can see, procedure Q1 converges but I couldn't get the resutls from Q2.

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

 

Sincerely yours

Amir

 

restart;

Eq1:=diff(f(x),x$3)+diff(f(x),x$2)*f(x)+b^2*sqrt(2*reynolds)*diff(diff(f(x),x$2)^2*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)^2*x,x$1);
eq1:=isolate(Eq1,diff(f(x),x,x,x));
eq2:=subs(g=f,isolate(Eq2,diff(g(x),x,x,x)));
EQ:=diff(f(x),x,x,x)=piecewise(x<c*0.1,rhs(eq1),rhs(eq2));
Eq11:=diff(theta(x),x$2)+pr*diff(theta(x),x$1)*f(x)+pr/prt*b^2*sqrt(2*reynolds)*diff(diff(f(x),x$2)*diff(theta(x),x$1)*x^2,x$1);
Eq22:=diff(g(x),x$2)+pr*diff(g(x),x$1)*f(x)+pr/prt*a^2*c*sqrt(2*reynolds)*diff(diff(f(x),x$2)*diff(g(x),x$1)*x^1,x$1);
eq11:=isolate(Eq11,diff(theta(x),x,x));
eq22:=subs(g=theta,isolate(Eq22,diff(g(x),x,x)));
EQT:=diff(theta(x),x,x)=piecewise(x<c*0.1,rhs(eq11),rhs(eq22));
EQT1a:=eval(EQT,EQ):
EQT2:=eval(EQT1a,{f(x)=G0(x),diff(f(x),x)=G1(x),diff(f(x),x,x)=G2(x)}):
bd:=c;
a:=0.13:
b:=0.41:
pr:=1;
prt:=0.86;
reynolds:=12734151.135786774055543653356602;     #10^6;   #1.125*10^8:

c:=88.419896050808975395120916434619:
;
Q:=proc(pp2) local res,F0,F1,F2;
print(pp2);
if not type(pp2,numeric) then return 'procname(_passed)' end if:
res:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=pp2},numeric,output=listprocedure);
F0,F1,F2:=op(subs(subs(res),[f(x),diff(f(x),x),diff(f(x),x,x)])):
F1(bd)-1;
end proc;
fsolve(Q(pp2)=0,pp2=(0..1002));
se:=%;
res2:=dsolve({EQ,f(0)=0,D(f)(0)=0,(D@@2)(f)(0)=se},numeric,output=listprocedure):
G0,G1,G2:=op(subs(subs(res2),[f(x),diff(f(x),x),diff(f(x),x,x)])):
plots:-odeplot(res2,[seq([x,diff(f(x),[x$i])],i=1..1)],0..c);



Q2:=proc(rr2) local solT,T0,T1;
print(rr2);
if not type(rr2,numeric) then return 'procname(_passed)' end if:
solT:=dsolve({EQT2,theta(0)=1,D(theta)(0)=-rr2},numeric,known=[G0,G1,G2],output=listprocedure):
T0,T1:=op(subs(subs(res),[theta(x),diff(theta(x),x)])):
T0(bd);
end proc;
fsolve(Q2(rr2)=0,rr2=(0..100));


shib:=%;
sol:=dsolve({EQT2,theta(0)=1,D(theta)(0)=-shib},numeric,known=[G0,G1,G2],output=listprocedure):
plots:-odeplot(sol,[x,theta(x)],0..c);
#fsolve(Q2(pp3)=0,pp3=-2..2):

Amir

Good day,

How can this be corrected ''Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging'' see the worksheet here VT.mw

Hi All,

I want to solve for G, but cant seem to figure it out. Can anyone tell me what am I missing or doing wrong?

 

SSD:= ((155/9)*2.5)+((155/9)^2/(2*9.807*0.346));
Lc:= 2*SSD-(2*(sqrt(1.08)+sqrt(0.6))^2)/(G+0.04);
Ls:= 2*SSD-(2*(0.6+SSD*tan(Pi/180)))/(G);
Lcon:= 202-Lc-Ls;
Eq:= (G-0.043)(Lc/2)+(-G)(Ls/2)+G*Lcon-0.043*Lc=14;


fsolve(Eq);

 

thank you

Hi all,

Thanks for helping me to solve the problem below using Maple.

dsys := {(1-4*(diff(ln(v(z)), z)))*(diff(u(z), z))+((3/2)*z^{-1}-2*(3* z^{-1} *(diff(ln(v(z)), z))+2*(diff(ln(v(z)),z,z )))))*u(z) = 0, -z*(diff(v(z), z))-v(z)+v(z)^(1/2)*u(z) = 0, v(0) = 1, u(0) = 1, (D(v))(0) = 1/4, (D(u))(0) = 3/8}

When trying    sol := dsolve(dsys, numeric)

I got : Error, (in DEtools/convertsys) unable to convert to an explicit first-order system.

 

Note that the analytic solution for z<=0 is:

if z>-4   then  u(z)=(1+z/4)exp(z/8) and   v(z)=exp(z/4)

else u(z)=0  and v(z)= (-4/z)exp(-1)

Regards


 \int_{a}^{b} f(x) \, dx \approx \tfrac{3h}{8}\left[f(x_0) + 3f(x_1) + 3f(x_2) 
+ 2f(x_3) + 3f(x_4) + 3f(x_5) + 2f(x_6) + ... + f(x_n)\right] .


A Fantastic Good morning Ladies and Gentelmen,

You have the mission to find the problem in that code/ give a better code?
I want to write a procedure for the simpson 3/8 rule using the above formula I took from wikipedia :
(I don't want to use any super maple command that estimate the simpson 3/8, but just using this formula and starting from scratch, I am a beginner so feel free to laugh at my code hehehe)
here is the code I managed to write after 365 hours :D

restart:
simp38:=proc(f,x0,xn,n)                 
     local h,summ,i,x,Integral:                                              
     h=evalf((xn-x0)/n);               #h is the distance between 2 points
     summ:=f(x0)+f(xn);              #initialise the summ
    
      for i from 1 to n-1 by 3 do             
          x[i]:=x0+i*h:
          summ:=evalf(summ+3*f(x[i])):
      od:
      for i from 2 to n-1 by 3 do              
          x[i]:=x0+i*h:
          summ:=evalf(summ+3*f(x[i])):
      od:
      for i from 3 to n-1 by 3 do                  
          x[i]:=x0+i*h:
          summ:=evalf(summ+2*f(x[i])):  
      od:
    
      Integral:=3*h/8*summ;                  
 end:

simp38(x->(x^5-5.15*x^2+8.55*x-4.05045),1,5,6);   #testing it      --->2481.087090
evalf(int(x^5-5.15*x^2+8.55*x-4.05045,x=1..5));      #comparing it --->2477.531533


Your help will be apreciated.


Hello!

So as the title says i have a problem with solution of second order differential equation system.

The problem is that after i set some initial conditions and check the solution when the time equals zero, the solution graph is completely different from my initial conditions which is obviously incorrect. I know that i can change "spacestep" parameter, but i am not sure if it is the only way to fix the problem. Are there any other ways to fix this?

Here is a part of my code if it helps

r0=1.570796327
theta0=6.000000000


Is it possible to use an option similar to range when using lsode method for dsolve? The ODEs I am trying to solve is stiff and will not work with the flag stiff=true or with method=rosenbrock unless i set Digits:=20. I want to avoid doing that as much as possible, since I believe wit will be very taxing, computationally. I have a very large systeom to solve. I found that method=lsode works with the default Digits=15. 

 

However I need to have the solutions in a given range stored for future access and manipulations. Using range gives me an error: 

Error, (in dsolve/numeric/an_args/lsode) lsode keyword was range, optional keyword must be one of 'ctrl', 'initial', 'itask', 'output', 'procedure', 'procvars', 'start', 'number', 'abserr', 'relerr', 'maxfun', 'minstep', 'maxstep', 'initstep', 'startinit', 'implicit', 'optimize', 'complex'

 

I cannot figure out how to use range or something similar with lsode. Anyone knows? 

Hi! I have to solve a very large system of ODEs for a set of functions that can be indexed with two integers, n and l, call them f[n,l](x), where n ranges from 1 to an order 1000 number and l ranges from 1 to an order 10 number. 

What is the best strategy in dsolve, keeping in mind I will only need to have quick access to the values of the first few such functions, say f[0,0](x) and f[0,1](x) for any given x in a given range. I dont need to store the rest, and trying to store all of them actually will give me a warning about length of output exceeds limit of 1000000. So, having maple output the result of dsolve in a procedural form seems to be what I want. However in doing so, for evey value of x i call the solving procedure maple will solve the system each time, thus taking a long time. 

 

My dilema is: to save storage space and memory I want to use dsolve output as procedures. To save time, when accessing the functions I need (only the first few in the large set of unknown functions that obey the ODE) for a given value of x, I would like maple to have those already stored instead of computing them with each invocation. So, what I need is a somehow split behavior of output. For some of my unknown functions to act as if i use range in dsolve, so that I have access at the values without further computation, and for the rest of the solutions (most of the unknown functions) just keep them in procedural form. How do I achieve this? 

 

Below is a schematic of my probem:

 

ODEs:={Set of about 10000 coupled ODEs for functions labeled f[n,l] n ranging from 1..1000 and l from 1..10}

 

dsolve(ODEs,numeric, output=???)#What options to put in dsolve such that I have quick access to f[0,0](x) and f[0,1](x) in a given range (say xmin=0.1 and xmax=10 for example) but not store the rest. 

 

 

Hello,

The idea: parameter "a" will have a new random value each 10 days.

The way I did it is working but it can get very long especially if I do it for a system of equations and for long time more than a year (365 days).

The code:

with(DEtools); with(plots);
n := 5;

for i to n do Ra[i] := RandomTools:-Generate(distribution(Uniform(0.1e-1, .5))); a[[i]] := Ra[i] end do;

b := 0.1e-2;

T := 10;

 eq := diff(L(t), t) = a*L(t)-b;

init[1] := L(0) = 100;
 sol[1] := dsolve({init[1], subs(a = a[[1]], eq)}, L(t), range = 0 .. T, numeric);


init[2] := L(T) = rhs(sol[1](T)[2]);

sol[2] := dsolve({init[2], subs(a = a[[2]], eq)}, L(t), range = T .. 2*T, numeric);

 

init[3] := L(2*T) = rhs(sol[2](2*T)[2]);
sol[3] := dsolve({init[3], subs(a = a[[3]], eq)}, L(t), range = 2*T .. 3*T, numeric);

p[1] := odeplot(sol[1], [t, L(t)], t = 0 .. T);

p[2] := odeplot(sol[2], [t, L(t)], t = T .. 2*T);

p[3] := odeplot(sol[3], [t, L(t)], t = 2*T .. 3*T);

p := display([p[1], p[2], p[3]]);
display(p);

 

Thank you

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