mehdi jafari

754 Reputation

13 Badges

11 years, 171 days

MaplePrimes Activity


These are replies submitted by mehdi jafari

@mmcdara here is the corrected code

abcdef.mw

@mmcdara 

Finger crossed! 
This code is amazing,I appreciate your efforts for solving this equation.
I have some offers to improve the code.
1- could you please put a “ minus “ within this prentice:
exp(nu*n*tau)
Is located in the following term:

G := cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*2*k[b]*T*w[c]^2
     *
     (w[c]*exp(-w[c]*tau)/(w[c]^2) + 2*Sum((w[c]*exp(-w[c]*tau)-nu*n*exp(nu*n*tau))/(w[c]^2-(nu*n)^2),n=1..5));

2- About the approximation, yes this is perfect.

3- This kind of differential integral equation in physics is about population of excited state in an atom, so it should be positive always, and it should fluctuate between 0.5 and 0, then I am thinking why this is negative and bigger than 0.5? What do you think?

@mmcdara When i run the your code, i got some errors in the code, is the code ok for you? If it is ok, then i explain the problem more.

As you're saying, this kind of differential integral equations are very complicated and i appreciate your efforts for solving this. I have some suggestions and explanations to improve this code,i hope we altogether with the help of you and other experts. could solve the problem.

1-# for n > 0 the denominator of op S is equivalent to -3.947841762*10^5*n^2 = (628.3185308*abs(n))^2 and the
# numerator is equivalent to 628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau):
# Thus, assuming n > 0:) but in this work for getting the following equation
  
k1:=2*k[b]*T*w[c]^2*sum((w[c]*exp(-w[c]*tau)-abs(nu[n])*exp(-abs(nu[n])*tau))/(w[c]^2-nu[n]^2),n=-infinity..infinity):

for calculating k1, we used a Taylor expansion which it has a sum on n=-infinity..infinity, and term of (nu[n]) is a frequency so it should be positive so we considered abs in k1 formula,  so I think we should consider n<0 as n>=0 ( n=-infinity..infinity, ). what do you think? how we can extend this code to consider entire domain?


2- the other approximation is considered in
;

sum((0.1000000000e-1*exp(-0.1000000000e-1*tau)-628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau))/(-394784.1762*n^2+0.1000000000e-3), n = -infinity .. infinity)

 I think so, it is true only for very small amounts in comparison the others, what do you think?is it general?

3-in equation (7)
RHS := eval(-G*rho(tau)+1/2(G-F), rho(tau)=U);
could you please put a * between 1/2 and (G-F)? and i think it can affect the output.


4-the initial value of rho=0.5;

5-the maximum value of t you want to solve to:
almost after 200 (or more)  it  'maybe'  be constant, for example for T=100 the final constant is 0.4975 (steady state).

thank you for the time you spend on this,
Sincerely yours

@mmcdara thank you, i think i should use numerical solution. tnx

@Carl Love can we remove  " t "  from expression between y and x? i mean can we find an explicit expression between y and x with removing the variable t ?

@tomleslie @Carl Love 
I mean when i assign stepsize to have accurate results, the results change significantly, such as:

 

restart:with(DEtools):

#constants
a0 := 10:sigma := 100: theta := (1/4)*Pi:L := 80000:omegace:=0:

#Differential Equations
Eq1:= diff(Rhoy(tau),tau)=-gam(tau)*(Z(tau)/(L))*Y(tau)+sin(theta)*(((1/2)*a0*exp(xi(tau)-400)/sigma^2)*diff(xi(tau),tau))-(omegace/omega)*Rhoz(tau):
Eq2:=diff(Rhoz(tau),tau)=(((1/2)*a0*exp(xi(tau)-400)/sigma^2)*cos(theta)^2+Rhoy(tau)*sin(theta))*((1/2)*a0*exp(xi(tau)-400)/sigma^2)*(1)+(omegace/omega)*Rhoy(tau):
Eq3:=diff(Y(tau),tau)=Rhoy(tau):
Eq4:=diff(Z(tau),tau)=Rhoz(tau):
Eq5:=diff(xi(tau),tau)=gam(tau)-Rhoz(tau):
Eq6:=diff(gam(tau), tau) = (1/4)*(2*a0^2*(exp(xi(tau)-400))^2*cos(theta)^2*(diff(xi(tau), tau))/sigma^4+8*Rhoy(tau)*(diff(Rhoy(tau), tau))+8*Rhoz(tau)*(diff(Rhoz(tau), tau)))/sqrt(4+a0^2*(exp(xi(tau)-400))^2*cos(theta)^2/sigma^4+4*Rhoy(tau)^2+4*Rhoz(tau)^2):

sys:=seq(Eq||i,i=1..6):

#initial conditions
inits:=xi(0)=400,gam(0)=1,Z(0)=0,Y(0)=0,Rhoz(0)=0,Rhoy(0)=0:

#sol:=dsolve([sys,inits]);

sol:=dsolve([sys,inits],numeric,method=dverk78,output=listprocedure):

GAM:=eval(gam(tau),sol):

XI:=eval(xi(tau),sol):

 

plot([GAM(tau)/51],tau=0..10,axes=boxed,gridlines,color=[blue,red]):
#GAM(tau) is between 0 to 1.2

 

sol1:=dsolve([sys,inits],numeric,stepsize=1e-5,output=listprocedure):#just add stepsize=1e-5

GAM:=eval(gam(tau),sol1):

XI:=eval(xi(tau),sol1):

plot([GAM(tau)/51],tau=0..10,axes=boxed,gridlines,color=[blue,red]);
#GAM(tau) is between 0 to 1.01960784314 to 1.01960784326 !!

 

 


 

Download ode_problem.mw

@tomleslie Thank you for the answer. Is t2 the minimum of t in the positive range which holds for inequality? If  not, how can i find the minimum of t in the positive range?

@Carl Love thank you for your answer and curiosity. I will share the correct conditions whenever possible. 

@Carl Love thank you carl for the time you spend on this. everytime i learn a lot form your knowlegde and personlaity.
Sincerely Yours.

@vv thank you for the clarification. now another problem arises.
Error, (in fsolve) procedures don't evaluate to numeric types. could please help me to solve the issue? sorry for the inconvenience.
Inflation-Inverse-1.mw

@vv i couldn't solve the problem could please check this code? tnx in advance
Inflation-Inverse-2.mw

@pik1432 
 

restart:with(LinearAlgebra):

printlevel:=1

1

(1)

x:=<x1,x2>:
W:=<<w11, w12>|<w21, w22>>:
Cal:= (W.x)^%T.(W.x):

eq1:=Array([seq](diff(Cal,t),t=x))^+;  # 1

Vector(2, {(1) = (2*w11*x1+2*w21*x2)*w11+(2*w12*x1+2*w22*x2)*w12, (2) = (2*w11*x1+2*w21*x2)*w21+(2*w12*x1+2*w22*x2)*w22})

(2)

eq2:=(2*W^%T.W.x);#2

Vector(2, {(1) = (2*w11^2+2*w12^2)*x1+(2*w11*w21+2*w12*w22)*x2, (2) = (2*w11*w21+2*w12*w22)*x1+(2*w21^2+2*w22^2)*x2})

(3)

seq(is(eq2[i]=eq1[i]),i=1..numelems(eq1))

true, true

(4)

 

 


 

Download 2.mw

@vv thank you

@Carl Love thanks for the comment. So the answer of the integration is undefined?


 

restart

f := ((1 - a)^2 + a^2*((1 - exp(-y))*(1 - exp(-x)) - 2 + exp(-x) + exp(-y)) + a*(2 - exp(-x) - exp(-y) + (1 - exp(-y))*(1 - exp(-x))))/(1 - a*exp(-x)*exp(-y))^3;

((1-a)^2+a^2*((1-exp(-y))*(1-exp(-x))-2+exp(-x)+exp(-y))+a*(2-exp(-x)-exp(-y)+(1-exp(-y))*(1-exp(-x))))/(1-a*exp(-x)*exp(-y))^3

(1)

a := 3/10:f:

int(f*exp(-x), x= 0..y + t, AllSolutions):

s := 2*(int((int(f*exp(-x)*exp(-y), x = 0 .. y + t,AllSolutions)), y = 0 .. infinity,AllSolutions));

s := 2*piecewise(Re(t) < 0, undefined, undefined)

(2)

 

 


 

Download stat2.mw

@Kitonum 
 

``

restart

f:=[evalf(solve(x^6-3*x-5))]

[1.451571465, .5973639664+1.275126159*I, -.7760033304+.9926461157*I, -1.094292737, -.7760033304-.9926461157*I, .5973639664-1.275126159*I]

(1)

remove(has,f,I);

[1.451571465, -1.094292737]

(2)

 

 

``


 

Download remove.mw

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