Advanced Search

Ola

30 Reputation

3 Badges

4 years, 142 days
Contact Ola

MaplePrimes Activity


These are questions asked by Ola

I want to integrate a function numerical...

Ola 30
parameters dsolve numeric
February 02 2025
1 2

Good day, 

I have this system of differential equations, which I solved numerically. 

Now I need to integrate the function y(t) numerically to get the corresponding expression for 

Z= exp(int y(t) dt ) (integrate from a given t to 0)

then I need to plot that system of differential equation vs (Z) with (t) as the parameter.

can you please help?

 eq1 := diff(x(t), t)-(1/6)*(6*x(t)^3*y(t)+(2*y(t)^2-2)*x(t)^2+3*y(t)*(z(t)-2)*x(t)-2*y(t)^2+2)*sqrt(3) = 0;

eq2 := diff(y(t), t)-(1/6)*(y(t)-1)*sqrt(3)*(y(t)+1)*(6*x(t)^2+2*y(t)*x(t)+3*z(t)-2) = 0;

eq3 := diff(z(t), t)-(1/3)*z(t)*sqrt(3)*(6*y(t)*x(t)^2+2*x(t)*y(t)^2+3*z(t)*y(t)-2*x(t)-3*y(t)) = 0;

eq4 := diff(s(t), t)-(1/3)*(y(t)-1)*sqrt(3)*(y(t)+1)*(6*x(t)^2+2*y(t)*x(t)+3*z(t)-2)*y(t) = 0;

eq5 := diff(Q(t), t)+sqrt(3)*((4/3-z(t)-2*x(t)^2)*y(t)+(2/3*(1-y(t)^2))*x(t))*Q(t)+2*((2/3)*y(t)+x(t))*P(t)/sqrt(3)-(2/3)*R(t) = 0;   
eq6 := diff(P(t), t)+(sqrt(3)*(1-z(t)-2*x(t)^2)*y(t)+2*(2-y(t)^2)*x(t)/sqrt(3))*P(t)+2*sqrt(3)*x(t)*Q(t)+(1/2)*R(t) = 0;       
eq7 := diff(R(t), t)-sqrt(3)*((3*x(t)^2+(3/2)*z(t)+4)*y(t)-(1/3*(1-3*y(t)^2))*x(t))*R(t)-z(t)*P(t) = 0;
eq8 := diff(T(t), t)-sqrt(3)*((6*x(t)^2+3*z(t)+1)*y(t)-(2*(1-y(t)^2))*x(t))*T(t)-(1-y(t)^2)*(Q(t)-(2/3)*P(t)) = 0;   
syst := {eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8};

ics := x(0) = .1, y(0) = .9, z(0) = .1, s(0) = .2, Q(0) = .1, P(0) = .9, R(0) = .1, T(0) = .2;
  

sol1 := dsolve({ics, op(syst)}, {P(t), Q(t), R(t), T(t), s(t), x(t), y(t), z(t)}, type = numeric, output = listprocedure);
 

 

Hello, I need help please, I need to use...

Ola 30 Maple
ode dsolve numeric
December 18 2024
2 5

I have this system

eq1 := diff(x(t), t)-(1/6)*(6*x(t)^3*y(t)+(2*y(t)^2-2)*x(t)^2+3*y(t)*(z(t)-2)*x(t)-2*y(t)^2+2)*sqrt(3) = 0;
                         
eq2 := diff(y(t), t)-(1/6)*(y(t)-1)*sqrt(3)*(y(t)+1)*(6*x(t)^2+2*y(t)*x(t)+3*z(t)-2) = 0;
                                    
eq3 := diff(z(t), t)-(1/3)*z(t)*sqrt(3)*(6*y(t)*x(t)^2+2*x(t)*y(t)^2+3*z(t)*y(t)-2*x(t)-3*y(t)) = 0;

I solved it numerically using these ics

ics := x(0) = -0.01, y(0) = .99, z(0) = 0.01

sol1 := dsolve({ics, op(syst)}, {x(t), y(t), z(t)}, type = numeric, output = listprocedure)

I need to use the x(t), y(t), z(t) as follows
  X :=  eval(x(t), sol1)
  Y :=  eval(y(t), sol1)

Z :=  eval(z(t), sol1)

to solve the following system for P(t), Q(t), R(t) numerically 
eq4 := diff(R(t), t)-P(t)*Z-(-2*(-Y^2+2)*X/sqrt(3)+sqrt(3)*(-2*X^2+Z+4/3)*Y)*R(t) = 0;
eq5 := diff(Q(t), t)-(2/3)*R(t)+2*((1/3)*Y+X)*P(t)/sqrt(3)-(-2*(-Y^2+2)*X/sqrt(3)+2*sqrt(3)*(X^2-(1/2)*Z-2/3)*X)*Q(t) = 0;
eq6 := diff(P(t), t)+(1/2)*R(t)+2*sqrt(3)*X*Q(t)+(2*(-Y^2+2)*X/sqrt(3)+sqrt(3)*(-2*X^2+Z+1)*Y)*P(t) = 0;

Any help please? 

Numeric integration of expressions of d...

Ola 30 Maple 2017
ode dsolve numeric
May 05 2021
1 0

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

 

How do I use numerical solutions of no...

Ola 30 Maple
ode homework dsolve
March 04 2021
1 6

I need to treat Eqs. (1.85-1.86) as a system of equations that should be solved simultaneously, and use numerical solutions in Eqs. (1.87-1.90)
Do you know if there are in-built commands that can do this?

Page 1 of 1