## 749 Reputation

11 years, 37 days

## the corrected code...

@mmcdara here is the corrected code

abcdef.mw

## thank you for such amazing effort, reall...

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?

## tnx for creative suggestions...

@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

## thank you...

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

## can we remove " t " from expression be...

@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 w...

@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 !!
 >

## minimum in the positive range...

@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?

## I will share if i find correct condition...

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

## tnx...

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

## thank you for the clearfication...

@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

## i couldn't solve the problem could pleas...

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

## is...

 > restart:with(LinearAlgebra):
 > printlevel:=1
 (1)
 > x:=: W:=<|>: Cal:= (W.x)^%T.(W.x):
 > eq1:=Array([seq](diff(Cal,t),t=x))^+;  # 1
 (2)
 > eq2:=(2*W^%T.W.x);#2
 (3)
 > seq(is(eq2[i]=eq1[i]),i=1..numelems(eq1))
 (4)
 >
 >

@vv thank you

## the answer of the integration is undefin...

@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 := 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));
 (2)
 >
 >

## another way...

 > restart
 > f:=[evalf(solve(x^6-3*x-5))]
 (1)
 > remove(has,f,I);
 (2)
 >
 >