Advanced Search

Question: Numeric integration of expressions of dependent variables derived from a system of nonlinear differential equations

Question:Numeric integration of expressions of dependent variables derived from a system of nonlinear differential equations

Ola 30 Maple 2017
ode dsolve numeric differential-equations
May 05 2021
2

Can anyone please help with this?  I don't know what is the problem in this code.

 

> 

restart; _EnvExplicit := true; with(DEtools); with(PDEtools); with(plots)

true

 (1)
> 

NULL
> 

Check*the*constraint
> 

eq1 := diff(X(z), z, z)+((diff(H(z), z))/H(z)-2/(1+z))*(diff(X(z), z))+exp(-sqrt(2/3)*X(z))*ZETA*(diff(Y(z), z))^2/(2*sqrt(6))-Y(z)*exp(-sqrt(2/3)*X(z))*(1-exp(-sqrt(2/3)*X(z)))/(2*H(z)^2*sqrt(6)*(1+z)^2) = 0
> 

eq2 := diff(Y(z), z, z)+((diff(H(z), z))/H(z)-2/(1+z))*(diff(Y(z), z))-(1/3)*sqrt(6)*(diff(X(z), z))*(diff(Y(z), z))+(1-(1/2)*exp(-sqrt(2/3)*X(z)))/(2*ZETA*H(z)^2*(1+z)^2) = 0
> 

eq3 := -2*H(z)*(1+z)*(diff(H(z), z))+3*H(z)^2+(1/2*((diff(X(z), z))^2+(1/2)*ZETA*exp(-sqrt(2/3)*X(z))*(diff(Y(z), z))^2))*H(z)^2*(1+z)^2-(1/4)*exp(-sqrt(2/3)*X(z))*Y(z)*(1-(1/2)*exp(-sqrt(2/3)*X(z))) = 0

-2*H(z)*(1+z)*(diff(H(z), z))+3*H(z)^2+((1/2)*(diff(X(z), z))^2+(1/4)*ZETA*exp(-(1/3)*6^(1/2)*X(z))*(diff(Y(z), z))^2)*H(z)^2*(1+z)^2-(1/4)*exp(-(1/3)*6^(1/2)*X(z))*Y(z)*(1-(1/2)*exp(-(1/3)*6^(1/2)*X(z))) = 0

 (2)
> 

NULL
> 

``
> 

NULL
> 

NULL
> 

diffux := solve(eq1, dif(H(z), z))
> 

diffur := solve(eq2, dif(H(z), z))
> 

NULL
> 

NULL
> 

NULL
> 

NULL
> 

NULL
> 

NULL
> 

``
> 

Numeric*code
> 

with(ListTools)
> 

Seed := randomize()
> 

with(stats[random])
> 

sys := {subs(ZETA = 10, eq1), subs(ZETA = 10, eq2), subs(ZETA = 10, eq3)}
> 

NULL
> 

vv := 10^100; savf := []; L := 0; Digits := 30; MM := []; MMM := []; N := 1; kkk := []; MM; []; MMM := []
> 

"for hh from 1 by 1 to 1 do:L:=0: F:=0: b:=1: t:=2:while (t>b)do:"
> 

n := 1.9; ZETA := 10; i := uniform[.60658, 1](N); ics := Y(i) = 1.8, (D(Y))(i) = 0, X(1) = .4, (D(X))(i) = 0, H(i) = .7; solution := dsolve({ics, op(subs(ZETA = 10, sys))}, {H(z), X(z), Y(z)}, type = numeric, method = bvp, output = listprocedure); m := [subs(solution, X(z))]; kk := [subs(solution, Y(z))]
> 

"if op(kk(0))>.7 then t:=(b)/(5): else t:=5*b :end if:end do: for T from i by 10^(-1) to 1 do ls[T] := sqrt(0.29995*(1+ T)^(3)+5*10^(-5)*(1+ T)^(4)+0.7): sr[T] := op(kk(T)): L := L+(sr[T]-ls[T])^2 end do:F:=sqrt(L): if F<vv then vv:=F; sol:=solution:MM:=[op(MM),[vv,op(kk(1))]]:hhh:=vv:Pa:=[g,i]:else vv:=vv: MMM:=[op(MMM),vv]:end if:end do:"

Error, (in unknown) solution is only defined in the range .88696222680754726..1.

  
> 

NULL
> 

``
> 

RH := subs(sol, Y(z)); Ry := subs(sol, X(z))

Error, invalid input: subs received sol, which is not valid for its 1st argument

  

Error, invalid input: subs received sol, which is not valid for its 1st argument

  
> 

1-Mrad(1)+Rx(1)+Ry(1)

1-Mrad(1)+Rx(1)+Ry(1)

 (3)
> 

with(plots)
> 

defHLCDM := sqrt(.29995*(1+z)^3+5*10^(-5)*(1+z)^4+.7)
> 

plotqLCDM := plot((1+z)*(diff(defHLCDM, v))/defHLCDM-1, z = 0 .. 10, color = green, legend = ["q LCDM"])
> 

plotqfR := [odeplot(sol, [z, Pa[1]*y(z)/(Pa[1]-1)+1], z = 0 .. 10, color = blue, legend = ["q R^n"])]; display(plotqLCDM, plotqfR)
> 

plotHLCDM := plot(sqrt(.29995*(1+z)^3+5*10^(-5)*(1+z)^4+.7), z = 0 .. Pa[2], color = green, legend = ["H LCDM"], axes = boxed)
> 

plotHfR := [odeplot(sol, [z, psi(z)], z = 0 .. Pa[2], color = blue, legend = ["H R^n"], axes = boxed)]; display(plotHLCDM, plotHfR)
 
> 

psi*curve*fit
> 

with(plots)
> 

with(ListTools)
> 

with(CurveFitting)
> 

points1 := PolynomialInterpolation([[0, RH(0)], [1, RH(1)], [10, RH(10)], [100, RH(100)], [400, RH(400)], [800, RH(800)], [1200, RH(1200)], [1600, RH(1600)], [2000, RH(2000)], [2400, RH(2400)], [2800, RH(2800)], [3200, RH(3200)], [3600, RH(3600)], [4000, RH(4000)], [4400, RH(4400)], [4600, RH(4600)], [4800, RH(4800)], [5000, RH(5000)], [5200, RH(5200)], [5400, RH(5400)], [5600, RH(5600)], [5800, RH(5800)], [6000, RH(6000)], [6200, RH(6200)], [Pa[2], RH(Pa[2])]], v)
 
> 

sd := plot(points1, z = 0 .. Pa[2], color = red, legend = ["curve fit"])

Error, (in plot) expecting a real constant as range endpoint but received Pa[2]

  
> 

points11 := [[0, RH(0)], [.1, RH(.1)], [1, RH(1)], [10, RH(10)], [100, RH(100)], [400, RH(400)], [800, RH(800)], [1200, RH(1200)], [1600, RH(1600)], [2000, RH(2000)], [2400, RH(2400)], [2800, RH(2800)], [3200, RH(3200)], [3600, RH(3600)], [4000, RH(4000)], [4400, RH(4400)], [4600, RH(4600)], [4800, RH(4800)], [5000, RH(5000)], [5200, RH(5200)], [5400, RH(5400)], [5600, RH(5600)], [5800, RH(5800)], [6000, RH(6000)], [6200, RH(6200)], [Pa[2], RH(Pa[2])]]
> 

sd := plot(points1, z = 0 .. Pa[2], color = red, legend = ["curve fit"], axes = boxed)
> 

pntplot1 := pointplot(points11, symbol = BOX)
> 

polycurve := PolynomialInterpolation(points11, z)
> 

polyplot := plot(polycurve, z = 0 .. Pa[2], axes = boxed)
> 

display([plotHLCDM, sd])

 
> 

``
> 

NULL
> 

``

 

 
> 

``
> 

Luminosity*distnace
> 

eq1h := diff(Hdi1(z), z)-1/points1 = 0
> 

eq2h := diff(Hdi2(z), z)-1/defHLCDM = 0
> 

sysh := {eq1h, eq2h}
> 

icsh := Hdi1(0) = 0, Hdi2(0) = 0; solutionh := dsolve({icsh, op(sysh)}, {Hdi1(z), Hdi2(z)}, type = numeric, output = listprocedure)
> 

plotLDistanceCFIT := [odeplot(solutionh, [z, (1+z)*Hdi1(z)], 0 .. Pa[2], color = blue, legend = ["Luminosity distance Curve fit"])]; plotLDistanceLCDM := [odeplot(solutionh, [z, (1+z)*Hdi2(z)], 0 .. Pa[2], color = red, legend = ["Luminosity distance LCDM"])]; display(plotLDistanceCFIT, plotLDistanceLCDM)

 
> 

``
> 

``
>   
> 

 

 

Download DistanceModulus.mw

 
Please Wait...