Advanced Search

tomleslie

13836 Reputation

20 Badges

14 years, 339 days
Contact tomleslie

MaplePrimes Activity


These are replies submitted by tomleslie

Calculations and some graphs...

tomleslie 13836
May 25 2023
0 0

@MaPal93 

for data sets 9..12 are in the attached

> 

restart:
> 

_local(gamma):

 #Define the assumptions ex-ante (variances as real and positive, correlations in between -1 and +1 and so on...) - or Maple wouldn't know

> 

assume(gamma::real, lambda__1::real, lambda__2::real, lambda__3::real, sigma__epsilon1::real, sigma__epsilon2::real, nu__0[1]::real, nu__0[2]::real, rho__u[1, 2]::real, rho__u[1, 3]::real, rho__u[2, 3]::real, rho__v[1, 2]::real, sigma__u[1]::real, sigma__u[2]::real, sigma__u[3]::real, sigma__v[1]::real, sigma__v[2]::real):
> 

assume(0 <= gamma, 0 <= lambda__1, 0 <= lambda__2, 0 <= lambda__3, 0 <= sigma__epsilon1, 0 <= sigma__epsilon2, 0 <= nu__0[1], 0 <= nu__0[2], -1 <= rho__u[1, 2] and rho__u[1, 2] <= 1, -1 <= rho__u[1, 3] and rho__u[1, 3] <= 1, -1 <= rho__u[2, 3] and rho__u[2, 3] <= 1, -1 <= rho__v[1, 2] and rho__v[1, 2] <= 1, 0 <= sigma__u[1], 0 <= sigma__u[2], 0 <= sigma__u[3], 0 <= sigma__v[1], 0 <= sigma__v[2]):


1. FIRST PART

 

> 

E := b*Exp_nu1_S + c*Exp_nu2_S + a:
> 

V := 2*Cov_u12*d*e + 2*Cov_u13*d*f + 2*Cov_u23*e*f + Var_u1*d^2 + Var_u2*e^2 + Var_u3*f^2 + b^2*Var_nu1_S + 2*b*c*Cov_nu12_S + c^2*Var_nu2_S:
> 

argmin := E - gamma/2*V:
> 

a := X__1*(-X__1*lambda__1 - nu__0[1]) + X__2*(-X__2*lambda__2 - nu__0[2]) + X__3*(-X__3*lambda__3 + (-nu__0[1] - nu__0[2])):
b := X__1 + X__3:
c := X__2 + X__3:
d := -X__1*lambda__1:
e := -X__2*lambda__2:
f := -X__3*lambda__3:

Exp_nu1_S := SIG1_m_nu01*theta__11 + SIG2_m_nu02*theta__12 + nu__0[1]:
Exp_nu2_S := SIG1_m_nu01*theta__21 + SIG2_m_nu02*theta__22 + nu__0[2]:
> 

FOC_1 := diff(argmin, X__1) = 0:
> 

FOC_2 := diff(argmin, X__2) = 0:
> 

FOC_3 := diff(argmin, X__3) = 0:
> 

strat := convert(solve({FOC_1, FOC_2, FOC_3}, {X__1, X__2, X__3}), list)
> 

# First strategy

Strategy1 := collect(rhs(strat[1]), [SIG1_m_nu01, SIG2_m_nu02]); length(Strategy1); paramsS1 := indets(Strategy1); type(Strategy1, linear(SIG1_m_nu01)); type(Strategy1, linear(SIG2_m_nu02)); `&alpha;__1` := simplify(eval(Strategy1, `~`[`=`]([SIG1_m_nu01, SIG2_m_nu02], 0)))
 
> 

# Second strategy


Strategy2 := collect(rhs(strat[2]), [SIG1_m_nu01, SIG2_m_nu02]); length(Strategy2); paramsS2 := indets(Strategy2); type(Strategy2, linear(SIG1_m_nu01)); type(Strategy2, linear(SIG2_m_nu02)); `&alpha;__2` := simplify(eval(Strategy2, `~`[`=`]([SIG1_m_nu01, SIG2_m_nu02], 0)))
> 

# Third strategy

Strategy3 := collect(rhs(strat[3]), [SIG1_m_nu01, SIG2_m_nu02]); length(Strategy3); paramsS3 := indets(Strategy3); type(Strategy3, linear(SIG1_m_nu01)); type(Strategy3, linear(SIG2_m_nu02)); `&alpha;__3` := simplify(eval(Strategy3, `~`[`=`]([SIG1_m_nu01, SIG2_m_nu02], 0)))


2. SECOND PART

 

> 

# Check linear projection theorem's calculations first:
eql1 := (beta__11*Var_nu1+beta__12*Cov_nu12)/(((beta__11)^2)*Var_S1+((beta__12)^2)*Var_S2+2*beta__11*beta__12*Cov_S12+Var_u1):

# Now plug in the betas:
beta__11 := simplify(coeff(Strategy1, SIG1_m_nu01)):
beta__12 := simplify(coeff(Strategy1, SIG2_m_nu02)):




# Check linear projection theorem's calculations first:
eql2 := (beta__21*Cov_nu12+beta__22*Var_nu2)/(((beta__21)^2)*Var_S1+((beta__22)^2)*Var_S2+2*beta__21*beta__22*Cov_S12+Var_u2):

# Now plug in the betas:
beta__21 := simplify(coeff(Strategy2, SIG1_m_nu01)):
beta__22 := simplify(coeff(Strategy2, SIG2_m_nu02)):

# Check linear projection theorem's calculations first:
eql3 := (beta__31*Var_nu1+(beta__31+beta__32)*Cov_nu12+beta__32*Var_nu2)/(((beta__31)^2)*Var_S1+((beta__32)^2)*Var_S2+2*beta__31*beta__32*Cov_S12+Var_u3):

# Now plug in the betas:
beta__31 := simplify(coeff(Strategy3, SIG1_m_nu01)):
beta__32 := simplify(coeff(Strategy3, SIG2_m_nu02)):
 
> 

MyEqs  := {eql1 - lambda__1, eql2 - lambda__2, eql3 - lambda__3}:
MyVars := {lambda__1, lambda__2, lambda__3}:

# [i = {....}] means: the equation number i contains this set of variables
print~([ seq([i=indets(MyEqs[i], name) intersect MyVars], i=1..nops(MyEqs)) ]):

[1 = {lambda__1, lambda__2, lambda__3}]

  

[2 = {lambda__1, lambda__2, lambda__3}]

  

[3 = {lambda__1, lambda__2, lambda__3}]

 (1)
> 

# General Calibration:

P := indets(MyEqs, name) minus MyVars:

d := 1/(((sigma__v[1])^2+(sigma__epsilon1)^2)*((sigma__v[2])^2+(sigma__epsilon2)^2)-(sigma__v[1])^2*(rho__v[1, 2])^2*(sigma__v[2])^2):

cal := [
  Cov_S12 = sigma__v[1]*rho__v[1, 2]*sigma__v[2],
  Cov_u12 = sigma__u[1]*rho__u[1, 2]*sigma__u[2],
  Cov_u13 = sigma__u[1]*rho__u[1, 3]*sigma__u[3],
  Cov_u23 = sigma__u[2]*rho__u[2, 3]*sigma__u[3],
  Var_S1 = (sigma__v[1])^2+(sigma__epsilon1)^2,
  Var_S2 = (sigma__v[2])^2+(sigma__epsilon2)^2,
  Var_nu1 = (sigma__v[1])^2,
  Var_nu2 = (sigma__v[2])^2,
  Var_u1 = (sigma__u[1])^2,
  Var_u2 = (sigma__u[2])^2,
  Var_u3 = (sigma__u[3])^2,
  Cov_nu12 = sigma__v[1]*rho__v[1, 2]*sigma__v[2],
  Cov_nu12_S = d*sigma__v[1]*rho__v[1, 2]*sigma__v[2]*(sigma__epsilon1)^2*(sigma__epsilon2)^2,
  Var_nu1_S = d*(sigma__epsilon1)^2*((sigma__v[1])^2*((sigma__v[2])^2+(sigma__epsilon2)^2)-(sigma__v[1])^2*(rho__v[1, 2])^2*(sigma__v[2])^2),
  Var_nu2_S = d*(sigma__epsilon2)^2*((sigma__v[2])^2*((sigma__v[1])^2+(sigma__epsilon1)^2)-(sigma__v[1])^2*(rho__v[1, 2])^2*(sigma__v[2])^2),
  theta__11 = 1-d*((sigma__epsilon1)^2*((sigma__v[2])^2+(sigma__epsilon2)^2)),
  theta__12 = d*sigma__v[1]*rho__v[1, 2]*sigma__v[2]*(sigma__epsilon1)^2,
  theta__21 = d*sigma__v[1]*rho__v[1, 2]*sigma__v[2]*(sigma__epsilon2)^2,
  theta__22 = 1-d*((sigma__epsilon2)^2*((sigma__v[1])^2+(sigma__epsilon1)^2))
]:

MyEqs_cal := eval(MyEqs,cal):
length(MyEqs_cal):

P_cal := indets(MyEqs_cal, name) minus MyVars:
> 

# Numerical Calibrations:

P_cal_cor := P_cal minus {gamma,sigma__epsilon1,sigma__epsilon2,sigma__v[1],sigma__v[2],sigma__u[1],sigma__u[2],sigma__u[3]}:

P_cal_var := P_cal minus {rho__v[1, 2],rho__u[1, 2],rho__u[1, 3],rho__u[2, 3]}:

zerocor := {rho__v[1, 2],rho__u[1, 2],rho__u[1, 3],rho__u[2, 3]}=~0:

normvar := {gamma,sigma__epsilon1,sigma__epsilon2,sigma__v[1],sigma__v[2],sigma__u[1],sigma__u[2],sigma__u[3]}=~1:

ncal1 := ({gamma = INDEX} union (P_cal_var minus {gamma} =~ 1)) union zerocor:
ncal2 := ({sigma__epsilon1 = INDEX} union (P_cal_var minus {sigma__epsilon1} =~ 1)) union zerocor:
ncal3 := ({sigma__epsilon2 = INDEX} union (P_cal_var minus {sigma__epsilon2} =~ 1)) union zerocor:
ncal4 := ({sigma__u[1] = INDEX} union (P_cal_var minus {sigma__u[1]} =~ 1)) union zerocor:
ncal5 := ({sigma__u[2] = INDEX} union (P_cal_var minus {sigma__u[2]} =~ 1)) union zerocor:
ncal6 := ({sigma__u[3] = INDEX} union (P_cal_var minus {sigma__u[3]} =~ 1)) union zerocor:
ncal7 := ({sigma__v[1] = INDEX} union (P_cal_var minus {sigma__v[1]} =~ 1)) union zerocor:
ncal8 := ({sigma__v[2] = INDEX} union (P_cal_var minus {sigma__v[2]} =~ 1)) union zerocor:
ncal9 := ({rho__u[1, 2] = INDEX} union (P_cal_cor minus {rho__u[1, 2]} =~ 0)) union normvar:
ncal10 := ({rho__u[1, 3] = INDEX} union (P_cal_cor minus {rho__u[1, 3]} =~ 0)) union normvar:
ncal11 := ({rho__u[2, 3] = INDEX} union (P_cal_cor minus {rho__u[2, 3]} =~ 0)) union normvar:
ncal12 := ({rho__v[1, 2] = INDEX} union (P_cal_cor minus {rho__v[1, 2]} =~ 0)) union normvar:
> 

  betseq:=["beta__11", "beta__12", "beta__21", "beta__22", "beta__31", "beta__32"]:
  lamseq:=["lambda1", "lambda2", "lambda3"]:
> 

#
# Data set 9
#
  arrInd:=1:
  lamRes9:=Array(1..1001, 1..4):
  betRes9:=Array(1..1001, 1..7):
  for INDEX from -1 by 0.002 to 1 do
      lam := fsolve( eval(MyEqs_cal, ncal9),
                     {lambda__1 = 0..1, lambda__2 = 0..1, lambda__3 = 0..1}
                   );
      betRes9[arrInd]:= Array( [ INDEX,
                                 seq(eval(eval(j, cal), lam union ncal9), j in parse~(betseq))
                               ]
                             ):
      lamRes9[arrInd]:= Array( [ INDEX,
                                 rhs~(convert(lam,list))[]
                               ]
                             ):
      arrInd:=arrInd+1;
  end do:
> 

#
# Data set 10
#
  arrInd:=1:
  lamRes10:=Array(1..1001, 1..4):
  betRes10:=Array(1..1001, 1..7):
  for INDEX from -1 by 0.002 to 1 do
      lam := fsolve( eval(MyEqs_cal, ncal10),
                     {lambda__1 = 0..1, lambda__2 = 0..1, lambda__3 = 0..1}
                   );
      betRes10[arrInd]:= Array( [ INDEX,
                                  seq(eval(eval(j, cal), lam union ncal10), j in parse~(betseq))
                                ]
                              ):
      lamRes10[arrInd]:= Array( [ INDEX,
                                  rhs~(convert(lam,list))[]
                                ]
                              ):
      arrInd:=arrInd+1;
  end do:
> 

#
# Data set 11
#
  arrInd:=1:
  lamRes11:=Array(1..1001, 1..4):
  betRes11:=Array(1..1001, 1..7):
  for INDEX from -1 by 0.002 to 1 do
      lam := fsolve( eval(MyEqs_cal, ncal11),
                     {lambda__1 = 0..1, lambda__2 = 0..1, lambda__3 = 0..1}
                   );
      betRes11[arrInd]:= Array( [ INDEX,
                                  seq(eval(eval(j, cal), lam union ncal11), j in parse~(betseq))
                                ]
                              ):
      lamRes11[arrInd]:= Array( [ INDEX,
                                  rhs~(convert(lam,list))[]
                                ]
                              ):
      arrInd:=arrInd+1;
  end do:
> 

#
# Data set 12
#
  arrInd:=1:
  lamRes12:=Array(1..1000, 1..4):
  betRes12:=Array(1..1000, 1..7):
  for INDEX from -0.998 by 0.002 to 1 do
      lam := fsolve( eval(MyEqs_cal, ncal12),
                     {lambda__1 = 0..1, lambda__2 = 0..1, lambda__3 = 0..1}
                   );
      betRes12[arrInd]:= Array( [ INDEX,
                                  seq(eval(eval(j, cal), lam union ncal12), j in parse~(betseq))
                                ]
                              ):
      lamRes12[arrInd]:= Array( [ INDEX,
                                  rhs~(convert(lam,list))[]
                                ]
                              ):
      arrInd:=arrInd+1;
  end do:
> 

  with(plots):
  cols:=[red, green, blue, cyan, black, orange]:
  lamplots:=Array(1..2, 1..2):
  betplots:=Array(1..2, 1..2):
  lamplots[1,1]:=display( seq( plot( lamRes9[..,1], lamRes9[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..4
                             ),
                          legend=lamseq,
                          title="lambda values - DataSet 9",
                          titlefont=[times, bold, 16]
                        ):
  lamplots[1,2]:=display( seq( plot( lamRes10[..,1], lamRes10[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..4
                             ),
                          legend=lamseq,
                          title="lambda values - DataSet 10",
                          titlefont=[times, bold, 16]
                        ):
  lamplots[2,1]:=display( seq( plot( lamRes11[..,1], lamRes11[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..4
                             ),
                          legend=lamseq,
                          title="lambda values - DataSet 11",
                          titlefont=[times, bold, 16]
                        ):
  lamplots[2,2]:=display( seq( plot( lamRes12[..,1], lamRes12[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..4
                             ),
                          legend=lamseq,
                          title="lambda values - DataSet 12",
                          titlefont=[times, bold, 16]
                        ):
  display( lamplots);
  betplots[1,1]:=display( seq( plot( betRes9[..,1], betRes9[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..7
                             ),
                          legend=betseq,
                          title="beta values - DataSet 9",
                          titlefont=[times, bold, 16]
                        ):
  betplots[1,2]:=display( seq( plot( betRes10[..,1], betRes10[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..7
                             ),
                          legend=betseq,
                          title="beta values - DataSet 10",
                          titlefont=[times, bold, 16]
                        ):
  betplots[2,1]:=display( seq( plot( betRes11[..,1], betRes11[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..7
                             ),
                          legend=betseq,
                          title="beta values - DataSet 11",
                          titlefont=[times, bold, 16]
                        ):
  betplots[2,2]:=display( seq( plot( betRes12[..,1], betRes12[..,j],
                                     color=cols[j-1]
                                   ),
                               j=2..7
                              ),
                          legend=betseq,
                          title="beta values - DataSet 12",
                          titlefont=[times, bold, 16]
                        ):
  display(betplots);

 
 

 

 

  

 

 
 

 

 

  

Download calPlots2.mw

You...

tomleslie 13836
May 19 2023
0 0

appear to have one equation and two dependent variables 'y' and 'u'. In general there should be the same number of equations and dependent variables.

Also the equation you did provide, appears to be first order in 'u', which means that you can only have one boundary condition on 'u'

You will need two boundary conditions on the dependent variable 'y'

You are doing it the hard way!!...

tomleslie 13836
May 19 2023
0 0

@The function 

For the second worksheet I wouldn't bother with plots:-inequal(), with or without polar coordinates. Check out the attached which produces the required plots using simple commands from the plottools() package.

> 

  restart;
  with(plots):
  with(plottools):
  display( sector( [0,0], 2, 0..Pi/2, color=red), scaling=constrained);
  display( sector( [0,0], 3, arcsin(sqrt(2)/2)..Pi/2+arcsin(sqrt(2)/2), color=red), scaling=constrained );

 

 
> 

 

Download simple2.mw

Ooooops...

tomleslie 13836
May 17 2023
0 0

I forgot that you wanted an "off-centre" rotation, which is actually a pretty trivial adjustment to the worksheet I posted previously. See the attached

> 

restart; with(plots); with(plottools); A := 2; p0 := display(sphere([0, 0, 0], .5, color = green, style = surface), seq(spacecurve([A*sin(phi)^2, j, phi], phi = 0 .. 2*Pi, coords = spherical, scaling = constrained, color = red, thickness = 4, linestyle = dot), j = 0 .. evalf(5*Pi*(1/6)), evalf((1/6)*Pi)), spacecurve({[sqrt(theta), 1, 1/10], [sqrt(theta), 2, 1/10], [sqrt(theta), 3, -1/10], [sqrt(theta), 1.7, 1/7]}, theta = 0 .. 2*Pi, coords = spherical, scaling = constrained, color = cyan, thickness = 4)); animate(rotate, [translate(p0, .3, 0, 0), 0, 0, j], j = 0 .. 2*Pi, frames = 20)

 
> 

 

> 

 

Download pul.mw

Upload a worksheet...

tomleslie 13836
May 17 2023
0 0

using the big green up-arrow in the Mapleprimes toolbar

A better(?) way...

tomleslie 13836
May 07 2023
0 0

to achieve the same outcome

> 

restart;
L:=[{{0, 1}, {0, 3}, {0, 5}, {0, 7}, {2, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 7}, {4, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 7}, {6, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 7}, {8, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {2, 7}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {4, 7}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {6, 7}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {7, 8}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {2, 5}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {4, 5}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {5, 6}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {5, 8}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {2, 3}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {3, 4}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {3, 6}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {3, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {3, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {5, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {7, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {8, 9}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {3, 6}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {5, 6}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {6, 7}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {6, 9}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {3, 4}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {4, 5}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {4, 7}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {4, 9}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 3}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 5}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 7}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 9}}, {{0, 1}, {2, 3}, {2, 5}, {2, 7}, {2, 9}}, {{0, 1}, {2, 3}, {3, 4}, {3, 6}, {3, 8}}, {{0, 1}, {2, 5}, {4, 5}, {5, 6}, {5, 8}}, {{0, 1}, {2, 7}, {4, 7}, {6, 7}, {7, 8}}, {{0, 1}, {2, 9}, {4, 9}, {6, 9}, {8, 9}}, {{0, 1}, {3, 4}, {4, 5}, {4, 7}, {4, 9}}, {{0, 1}, {3, 6}, {5, 6}, {6, 7}, {6, 9}}, {{0, 1}, {3, 8}, {5, 8}, {7, 8}, {8, 9}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 2}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 4}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 6}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 8}}, {{0, 3}, {1, 2}, {1, 4}, {1, 6}, {1, 8}}, {{0, 3}, {1, 2}, {2, 5}, {2, 7}, {2, 9}}, {{0, 3}, {1, 2}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {1, 4}, {2, 3}, {3, 6}, {3, 8}}, {{0, 3}, {1, 4}, {4, 5}, {4, 7}, {4, 9}}, {{0, 3}, {1, 6}, {2, 3}, {3, 4}, {3, 8}}, {{0, 3}, {1, 6}, {5, 6}, {6, 7}, {6, 9}}, {{0, 3}, {1, 8}, {2, 3}, {3, 4}, {3, 6}}, {{0, 3}, {1, 8}, {5, 8}, {7, 8}, {8, 9}}, {{0, 3}, {2, 3}, {3, 4}, {3, 6}, {5, 8}}, {{0, 3}, {2, 3}, {3, 4}, {3, 6}, {7, 8}}, {{0, 3}, {2, 3}, {3, 4}, {3, 6}, {8, 9}}, {{0, 3}, {2, 3}, {3, 4}, {3, 8}, {5, 6}}, {{0, 3}, {2, 3}, {3, 4}, {3, 8}, {6, 7}}, {{0, 3}, {2, 3}, {3, 4}, {3, 8}, {6, 9}}, {{0, 3}, {2, 3}, {3, 6}, {3, 8}, {4, 5}}, {{0, 3}, {2, 3}, {3, 6}, {3, 8}, {4, 7}}, {{0, 3}, {2, 3}, {3, 6}, {3, 8}, {4, 9}}, {{0, 3}, {2, 5}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {2, 5}, {4, 5}, {5, 6}, {5, 8}}, {{0, 3}, {2, 7}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {2, 7}, {4, 7}, {6, 7}, {7, 8}}, {{0, 3}, {2, 9}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {2, 9}, {4, 9}, {6, 9}, {8, 9}}, {{0, 5}, {1, 2}, {1, 4}, {1, 6}, {1, 8}}, {{0, 5}, {1, 2}, {2, 3}, {2, 7}, {2, 9}}, {{0, 5}, {1, 2}, {4, 5}, {5, 6}, {5, 8}}, {{0, 5}, {1, 4}, {2, 5}, {5, 6}, {5, 8}}, {{0, 5}, {1, 4}, {3, 4}, {4, 7}, {4, 9}}, {{0, 5}, {1, 6}, {2, 5}, {4, 5}, {5, 8}}, {{0, 5}, {1, 6}, {3, 6}, {6, 7}, {6, 9}}, {{0, 5}, {1, 8}, {2, 5}, {4, 5}, {5, 6}}, {{0, 5}, {1, 8}, {3, 8}, {7, 8}, {8, 9}}, {{0, 5}, {2, 3}, {3, 4}, {3, 6}, {3, 8}}, {{0, 5}, {2, 3}, {4, 5}, {5, 6}, {5, 8}}, {{0, 5}, {2, 5}, {3, 4}, {5, 6}, {5, 8}}, {{0, 5}, {2, 5}, {3, 6}, {4, 5}, {5, 8}}, {{0, 5}, {2, 5}, {3, 8}, {4, 5}, {5, 6}}, {{0, 5}, {2, 5}, {4, 5}, {5, 6}, {7, 8}}, {{0, 5}, {2, 5}, {4, 5}, {5, 6}, {8, 9}}, {{0, 5}, {2, 5}, {4, 5}, {5, 8}, {6, 7}}, {{0, 5}, {2, 5}, {4, 5}, {5, 8}, {6, 9}}, {{0, 5}, {2, 5}, {4, 7}, {5, 6}, {5, 8}}, {{0, 5}, {2, 5}, {4, 9}, {5, 6}, {5, 8}}, {{0, 5}, {2, 7}, {4, 5}, {5, 6}, {5, 8}}, {{0, 5}, {2, 7}, {4, 7}, {6, 7}, {7, 8}}, {{0, 5}, {2, 9}, {4, 5}, {5, 6}, {5, 8}}, {{0, 5}, {2, 9}, {4, 9}, {6, 9}, {8, 9}}, {{0, 7}, {1, 2}, {1, 4}, {1, 6}, {1, 8}}, {{0, 7}, {1, 2}, {2, 3}, {2, 5}, {2, 9}}, {{0, 7}, {1, 2}, {4, 7}, {6, 7}, {7, 8}}, {{0, 7}, {1, 4}, {2, 7}, {6, 7}, {7, 8}}, {{0, 7}, {1, 4}, {3, 4}, {4, 5}, {4, 9}}, {{0, 7}, {1, 6}, {2, 7}, {4, 7}, {7, 8}}, {{0, 7}, {1, 6}, {3, 6}, {5, 6}, {6, 9}}, {{0, 7}, {1, 8}, {2, 7}, {4, 7}, {6, 7}}, {{0, 7}, {1, 8}, {3, 8}, {5, 8}, {8, 9}}, {{0, 7}, {2, 3}, {3, 4}, {3, 6}, {3, 8}}, {{0, 7}, {2, 3}, {4, 7}, {6, 7}, {7, 8}}, {{0, 7}, {2, 5}, {4, 5}, {5, 6}, {5, 8}}, {{0, 7}, {2, 5}, {4, 7}, {6, 7}, {7, 8}}, {{0, 7}, {2, 7}, {3, 4}, {6, 7}, {7, 8}}, {{0, 7}, {2, 7}, {3, 6}, {4, 7}, {7, 8}}, {{0, 7}, {2, 7}, {3, 8}, {4, 7}, {6, 7}}, {{0, 7}, {2, 7}, {4, 5}, {6, 7}, {7, 8}}, {{0, 7}, {2, 7}, {4, 7}, {5, 6}, {7, 8}}, {{0, 7}, {2, 7}, {4, 7}, {5, 8}, {6, 7}}, {{0, 7}, {2, 7}, {4, 7}, {6, 7}, {8, 9}}, {{0, 7}, {2, 7}, {4, 7}, {6, 9}, {7, 8}}, {{0, 7}, {2, 7}, {4, 9}, {6, 7}, {7, 8}}, {{0, 7}, {2, 9}, {4, 7}, {6, 7}, {7, 8}}, {{0, 7}, {2, 9}, {4, 9}, {6, 9}, {8, 9}}, {{0, 9}, {1, 2}, {1, 4}, {1, 6}, {1, 8}}, {{0, 9}, {1, 2}, {2, 3}, {2, 5}, {2, 7}}, {{0, 9}, {1, 2}, {4, 9}, {6, 9}, {8, 9}}, {{0, 9}, {1, 4}, {2, 9}, {6, 9}, {8, 9}}, {{0, 9}, {1, 4}, {3, 4}, {4, 5}, {4, 7}}, {{0, 9}, {1, 6}, {2, 9}, {4, 9}, {8, 9}}, {{0, 9}, {1, 6}, {3, 6}, {5, 6}, {6, 7}}, {{0, 9}, {1, 8}, {2, 9}, {4, 9}, {6, 9}}, {{0, 9}, {1, 8}, {3, 8}, {5, 8}, {7, 8}}, {{0, 9}, {2, 3}, {3, 4}, {3, 6}, {3, 8}}, {{0, 9}, {2, 3}, {4, 9}, {6, 9}, {8, 9}}, {{0, 9}, {2, 5}, {4, 5}, {5, 6}, {5, 8}}, {{0, 9}, {2, 5}, {4, 9}, {6, 9}, {8, 9}}, {{0, 9}, {2, 7}, {4, 7}, {6, 7}, {7, 8}}, {{0, 9}, {2, 7}, {4, 9}, {6, 9}, {8, 9}}, {{0, 9}, {2, 9}, {3, 4}, {6, 9}, {8, 9}}, {{0, 9}, {2, 9}, {3, 6}, {4, 9}, {8, 9}}, {{0, 9}, {2, 9}, {3, 8}, {4, 9}, {6, 9}}, {{0, 9}, {2, 9}, {4, 5}, {6, 9}, {8, 9}}, {{0, 9}, {2, 9}, {4, 7}, {6, 9}, {8, 9}}, {{0, 9}, {2, 9}, {4, 9}, {5, 6}, {8, 9}}, {{0, 9}, {2, 9}, {4, 9}, {5, 8}, {6, 9}}, {{0, 9}, {2, 9}, {4, 9}, {6, 7}, {8, 9}}, {{0, 9}, {2, 9}, {4, 9}, {6, 9}, {7, 8}}, {{1, 2}, {2, 3}, {2, 5}, {2, 7}, {4, 9}}, {{1, 2}, {2, 3}, {2, 5}, {2, 7}, {6, 9}}, {{1, 2}, {2, 3}, {2, 5}, {2, 7}, {8, 9}}, {{1, 2}, {2, 3}, {2, 5}, {2, 9}, {4, 7}}, {{1, 2}, {2, 3}, {2, 5}, {2, 9}, {6, 7}}, {{1, 2}, {2, 3}, {2, 5}, {2, 9}, {7, 8}}, {{1, 2}, {2, 3}, {2, 7}, {2, 9}, {4, 5}}, {{1, 2}, {2, 3}, {2, 7}, {2, 9}, {5, 6}}, {{1, 2}, {2, 3}, {2, 7}, {2, 9}, {5, 8}}, {{1, 2}, {2, 5}, {2, 7}, {2, 9}, {3, 4}}, {{1, 2}, {2, 5}, {2, 7}, {2, 9}, {3, 6}}, {{1, 2}, {2, 5}, {2, 7}, {2, 9}, {3, 8}}, {{1, 2}, {3, 4}, {4, 5}, {4, 7}, {4, 9}}, {{1, 2}, {3, 6}, {5, 6}, {6, 7}, {6, 9}}, {{1, 2}, {3, 8}, {5, 8}, {7, 8}, {8, 9}}, {{1, 4}, {2, 3}, {2, 5}, {2, 7}, {2, 9}}, {{1, 4}, {2, 3}, {4, 5}, {4, 7}, {4, 9}}, {{1, 4}, {2, 5}, {3, 4}, {4, 7}, {4, 9}}, {{1, 4}, {2, 7}, {3, 4}, {4, 5}, {4, 9}}, {{1, 4}, {2, 9}, {3, 4}, {4, 5}, {4, 7}}, {{1, 4}, {3, 4}, {4, 5}, {4, 7}, {6, 9}}, {{1, 4}, {3, 4}, {4, 5}, {4, 7}, {8, 9}}, {{1, 4}, {3, 4}, {4, 5}, {4, 9}, {6, 7}}, {{1, 4}, {3, 4}, {4, 5}, {4, 9}, {7, 8}}, {{1, 4}, {3, 4}, {4, 7}, {4, 9}, {5, 6}}, {{1, 4}, {3, 4}, {4, 7}, {4, 9}, {5, 8}}, {{1, 4}, {3, 6}, {4, 5}, {4, 7}, {4, 9}}, {{1, 4}, {3, 6}, {5, 6}, {6, 7}, {6, 9}}, {{1, 4}, {3, 8}, {4, 5}, {4, 7}, {4, 9}}, {{1, 4}, {3, 8}, {5, 8}, {7, 8}, {8, 9}}, {{1, 6}, {2, 3}, {2, 5}, {2, 7}, {2, 9}}, {{1, 6}, {2, 3}, {5, 6}, {6, 7}, {6, 9}}, {{1, 6}, {2, 5}, {3, 6}, {6, 7}, {6, 9}}, {{1, 6}, {2, 7}, {3, 6}, {5, 6}, {6, 9}}, {{1, 6}, {2, 9}, {3, 6}, {5, 6}, {6, 7}}, {{1, 6}, {3, 4}, {4, 5}, {4, 7}, {4, 9}}, {{1, 6}, {3, 4}, {5, 6}, {6, 7}, {6, 9}}, {{1, 6}, {3, 6}, {4, 5}, {6, 7}, {6, 9}}, {{1, 6}, {3, 6}, {4, 7}, {5, 6}, {6, 9}}, {{1, 6}, {3, 6}, {4, 9}, {5, 6}, {6, 7}}, {{1, 6}, {3, 6}, {5, 6}, {6, 7}, {8, 9}}, {{1, 6}, {3, 6}, {5, 6}, {6, 9}, {7, 8}}, {{1, 6}, {3, 6}, {5, 8}, {6, 7}, {6, 9}}, {{1, 6}, {3, 8}, {5, 6}, {6, 7}, {6, 9}}, {{1, 6}, {3, 8}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {2, 3}, {2, 5}, {2, 7}, {2, 9}}, {{1, 8}, {2, 3}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {2, 5}, {3, 8}, {7, 8}, {8, 9}}, {{1, 8}, {2, 7}, {3, 8}, {5, 8}, {8, 9}}, {{1, 8}, {2, 9}, {3, 8}, {5, 8}, {7, 8}}, {{1, 8}, {3, 4}, {4, 5}, {4, 7}, {4, 9}}, {{1, 8}, {3, 4}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {3, 6}, {5, 6}, {6, 7}, {6, 9}}, {{1, 8}, {3, 6}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {3, 8}, {4, 5}, {7, 8}, {8, 9}}, {{1, 8}, {3, 8}, {4, 7}, {5, 8}, {8, 9}}, {{1, 8}, {3, 8}, {4, 9}, {5, 8}, {7, 8}}, {{1, 8}, {3, 8}, {5, 6}, {7, 8}, {8, 9}}, {{1, 8}, {3, 8}, {5, 8}, {6, 7}, {8, 9}}, {{1, 8}, {3, 8}, {5, 8}, {6, 9}, {7, 8}}]:

  L1:=ListTools:-Categorize((x,y) -> x[1] = y[1], L):
#
# Check a few entries in the output list
#
  L1[1];
  L1[2];
  L1[9];

[{{0, 1}, {0, 3}, {0, 5}, {0, 7}, {2, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 7}, {4, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 7}, {6, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 7}, {8, 9}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {2, 7}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {4, 7}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {6, 7}}, {{0, 1}, {0, 3}, {0, 5}, {0, 9}, {7, 8}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {2, 5}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {4, 5}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {5, 6}}, {{0, 1}, {0, 3}, {0, 7}, {0, 9}, {5, 8}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {2, 3}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {3, 4}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {3, 6}}, {{0, 1}, {0, 5}, {0, 7}, {0, 9}, {3, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {3, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {5, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {7, 8}}, {{0, 1}, {1, 2}, {1, 4}, {1, 6}, {8, 9}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {3, 6}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {5, 6}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {6, 7}}, {{0, 1}, {1, 2}, {1, 4}, {1, 8}, {6, 9}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {3, 4}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {4, 5}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {4, 7}}, {{0, 1}, {1, 2}, {1, 6}, {1, 8}, {4, 9}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 3}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 5}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 7}}, {{0, 1}, {1, 4}, {1, 6}, {1, 8}, {2, 9}}, {{0, 1}, {2, 3}, {2, 5}, {2, 7}, {2, 9}}, {{0, 1}, {2, 3}, {3, 4}, {3, 6}, {3, 8}}, {{0, 1}, {2, 5}, {4, 5}, {5, 6}, {5, 8}}, {{0, 1}, {2, 7}, {4, 7}, {6, 7}, {7, 8}}, {{0, 1}, {2, 9}, {4, 9}, {6, 9}, {8, 9}}, {{0, 1}, {3, 4}, {4, 5}, {4, 7}, {4, 9}}, {{0, 1}, {3, 6}, {5, 6}, {6, 7}, {6, 9}}, {{0, 1}, {3, 8}, {5, 8}, {7, 8}, {8, 9}}]

  

[{{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 2}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 4}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 6}}, {{0, 3}, {0, 5}, {0, 7}, {0, 9}, {1, 8}}, {{0, 3}, {1, 2}, {1, 4}, {1, 6}, {1, 8}}, {{0, 3}, {1, 2}, {2, 5}, {2, 7}, {2, 9}}, {{0, 3}, {1, 2}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {1, 4}, {2, 3}, {3, 6}, {3, 8}}, {{0, 3}, {1, 4}, {4, 5}, {4, 7}, {4, 9}}, {{0, 3}, {1, 6}, {2, 3}, {3, 4}, {3, 8}}, {{0, 3}, {1, 6}, {5, 6}, {6, 7}, {6, 9}}, {{0, 3}, {1, 8}, {2, 3}, {3, 4}, {3, 6}}, {{0, 3}, {1, 8}, {5, 8}, {7, 8}, {8, 9}}, {{0, 3}, {2, 3}, {3, 4}, {3, 6}, {5, 8}}, {{0, 3}, {2, 3}, {3, 4}, {3, 6}, {7, 8}}, {{0, 3}, {2, 3}, {3, 4}, {3, 6}, {8, 9}}, {{0, 3}, {2, 3}, {3, 4}, {3, 8}, {5, 6}}, {{0, 3}, {2, 3}, {3, 4}, {3, 8}, {6, 7}}, {{0, 3}, {2, 3}, {3, 4}, {3, 8}, {6, 9}}, {{0, 3}, {2, 3}, {3, 6}, {3, 8}, {4, 5}}, {{0, 3}, {2, 3}, {3, 6}, {3, 8}, {4, 7}}, {{0, 3}, {2, 3}, {3, 6}, {3, 8}, {4, 9}}, {{0, 3}, {2, 5}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {2, 5}, {4, 5}, {5, 6}, {5, 8}}, {{0, 3}, {2, 7}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {2, 7}, {4, 7}, {6, 7}, {7, 8}}, {{0, 3}, {2, 9}, {3, 4}, {3, 6}, {3, 8}}, {{0, 3}, {2, 9}, {4, 9}, {6, 9}, {8, 9}}]

  

[{{1, 8}, {2, 3}, {2, 5}, {2, 7}, {2, 9}}, {{1, 8}, {2, 3}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {2, 5}, {3, 8}, {7, 8}, {8, 9}}, {{1, 8}, {2, 7}, {3, 8}, {5, 8}, {8, 9}}, {{1, 8}, {2, 9}, {3, 8}, {5, 8}, {7, 8}}, {{1, 8}, {3, 4}, {4, 5}, {4, 7}, {4, 9}}, {{1, 8}, {3, 4}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {3, 6}, {5, 6}, {6, 7}, {6, 9}}, {{1, 8}, {3, 6}, {5, 8}, {7, 8}, {8, 9}}, {{1, 8}, {3, 8}, {4, 5}, {7, 8}, {8, 9}}, {{1, 8}, {3, 8}, {4, 7}, {5, 8}, {8, 9}}, {{1, 8}, {3, 8}, {4, 9}, {5, 8}, {7, 8}}, {{1, 8}, {3, 8}, {5, 6}, {7, 8}, {8, 9}}, {{1, 8}, {3, 8}, {5, 8}, {6, 7}, {8, 9}}, {{1, 8}, {3, 8}, {5, 8}, {6, 9}, {7, 8}}]

 (1)
> 

 

Download listlist2.mw

Before...

tomleslie 13836
May 05 2023
0 0

@KIRAN SAJJAN 

you start worrying about plotting something, maybe you should take a look at the values yoo are generating!!

  1. All values of 'Theta1', apart from 'Theta1[11]' are of order 10-40. Do you think this is correct?
  2. All values of 'U1', apart from 'U1[11]' are of order 10-38. Do you think this is correct?

I should probably also point out...

tomleslie 13836
May 04 2023
0 0

that after a quick experiment Maple's built-in FDM method, cannot meet its own error criterion with the dsolve() maxmesh option set to 32 or 64, so tehe default setiing of 128 seems a pretty good choice for this problem although obviously you can increase it (all the way to 134217728, I believe).

Good luck with a 10-point mesh!!

Still trying - huh!...

tomleslie 13836
May 04 2023
0 0

@KIRAN SAJJAN 

because you don't like the (default) 128-point FDM method built-in to Maple? You think you can do better writinng your own 10-point FDM method?

Anyway, here are a few things you are going to have to fix - there are probably others

  1. One reason that this isn't even going to execute properly is that your definition of 'ode4' contains references to quantities such as 'Theta1(Y)' 'diff(U0(Y), Y)', 'diff(Theta0(Y), Y, Y)' etc, when these quantities should be replaced by the appropriate approximations given in d1, d2,.... d8.
  2. Another reason this code block isn't going to work is that values for 'R' are undefined.
  3. Another reason is that I am pretty sure you are going to have to define values for U0[-1], U1[-1], Theta0[-1] and Theta1[-1]

 

As I said in...

tomleslie 13836
April 26 2023
0 0

my previous post (which you should read!!)

I also had to use the DirectSearch package from the Maple Application centre in order to solve the reulting algebraic equations. You might be able to get a simple fsolve() to work with these equations by restricting the solution intervals for all variables. I don't have the time/patience to try this

In the attached I have generated values for f(eta) and Theta(eta) over the range 0..10 in steps of 0.5. You can compare these with the values in the variables data1 and data2. The agreement won't be very good for several reasons

  1. The Maple dsolve() command is generating a solution over the range 0..10
  2. Your FDM implementation is generating a solution over the range 0..5 (but using boundary conditions for eta=10??). You might get better agreement by changing your step size 'h' to 1.
  3. By default Maple's dsolve command uses a 128-point mesh (which will automatically increase if certain error criteria are not met). Your FDM implementation is using a 10-point mesh

NULL

> 

restart;
with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

 (1)
> 

A := 0.5;
M := 0.5;
R := 0.5;
Pr := 0.72;
Ec := 0.1;
g := 0.1;

.5

  

.5

  

.5

  

.72

  

.1

  

.1

 (2)
> 

ODE := diff(f(eta), eta, eta, eta) + f(eta)*diff(f(eta), eta, eta) - diff(f(eta), eta)^2 - M*diff(f(eta), eta) - A*diff(f(eta), eta) - A*eta*diff(f(eta), eta, eta)/2 = 0, (1 + R)*diff(Theta(eta), eta, eta) + Pr*(f(eta)*diff(Theta(eta), eta) - Theta(eta)*diff(f(eta), eta)) - Pr*A*(Theta(eta) + eta/2*diff(Theta(eta), eta)) + Pr*(Ec*diff(f(eta), eta, eta)^2 + g*Theta(eta));
BCS:=f(0)=0,D(f)(0)=1,Theta(0)=1,D(f)(10)=0,D(Theta)(10)=0;
sol:=dsolve({ODE,BCS},numeric,output = listprocedure):
p1:=odeplot(sol,[eta,f(eta)],0..9,color=green,linestyle=3,thickness=3,legend = ["dsolve"]):
p2:=odeplot(sol, [eta, Theta(eta)], 0 .. 9, color = green, linestyle = 3, thickness = 3, legend = ["dsolve"]):

diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(diff(f(eta), eta))^2-1.0*(diff(f(eta), eta))-.2500000000*eta*(diff(diff(f(eta), eta), eta)) = 0, 1.5*(diff(diff(Theta(eta), eta), eta))+.72*f(eta)*(diff(Theta(eta), eta))-.72*Theta(eta)*(diff(f(eta), eta))-.288*Theta(eta)-.1800000000*eta*(diff(Theta(eta), eta))+0.72e-1*(diff(diff(f(eta), eta), eta))^2

  

f(0) = 0, (D(f))(0) = 1, Theta(0) = 1, (D(f))(10) = 0, (D(Theta))(10) = 0

 (3)
> 

 fVals:= Matrix([[eta, f(eta)], seq( [j, eval(f(eta), sol)(j)], j=0..10, 0.5)]);
 ThetaVals:=Matrix([[eta, Theta(eta)], seq( [j, eval(Theta(eta), sol)(j)], j=0..10, 0.5)]);

_rtable[36893488148154145476]

  

_rtable[36893488148154135244]

 (4)
> 

d1:= (f[i+1]-f[i-1])/(2*h):
d2:=(f[i+1]-2*f[i]+f[i-1])/h^2:
d3:= (f[i+1]-3*f[i]+3*f[i-1]-f[i-2])/h^3:
d4:= (Theta[i+1]-Theta[i-1])/(2*h):
d5:=(Theta[i+1]-2*Theta[i]+Theta[i-1])/h^2:
h := 0.5:
ode1 := d3 + f[i]*d2 - d1^2 - M*d1 - A*d1 - A*eta*d2/2 = 0:
ode2 := (1 + R)*d5 + Pr*(-d1*Theta[i] + d4*f[i]) - Pr*A*(Theta[i] + eta/2*d4) + Pr*(Ec*d2^2 + g*Theta[i]):
> 

eta[0]:=0:
f[0]:=0:
Theta[0] := 1:
f[-1]:=f[1]-2*h:
f[10]:=f[9]:
Theta[10] := Theta[9]:
for i from 0 to 9 do
    eta[i+1]:=eta[i]+h;
    Eq[i]:=eval(ode1,eta=eta[i]);
    eq[i]:=eval(ode2,eta=eta[i]);
end do:
> 

 soln1:= DirectSearch:-SolveEquations( [Eq[0],Eq[1],Eq[2],Eq[3],Eq[4],Eq[5],Eq[6],Eq[7],Eq[8],Eq[9]])[3];
 soln2:= DirectSearch:-SolveEquations( eval([eq[0],eq[1],eq[2],eq[3],eq[4],eq[5],eq[6],eq[7],eq[8],eq[9]], soln1))[3];
 data1:=[seq([eta[0]+i*h,eval(f[i], soln1 )],i=0..9)];
 data2 := [seq([h*i + eta[0], eval(Theta[i], soln2)], i = 0 .. 9)];

[f[-2] = HFloat(-1.8128916235442254), f[1] = HFloat(0.3592770908855222), f[2] = HFloat(0.5657329830840051), f[3] = HFloat(0.6872174854707988), f[4] = HFloat(0.76034019145543), f[5] = HFloat(0.8052695951373906), f[6] = HFloat(0.8333150468603043), f[7] = HFloat(0.8508739043275494), f[8] = HFloat(0.8614824212622565), f[9] = HFloat(0.8667958074968981)]

  

[Theta[-1] = HFloat(1.425809612677366), Theta[1] = HFloat(0.7269817315917098), Theta[2] = HFloat(0.5464778444851895), Theta[3] = HFloat(0.4238106210804488), Theta[4] = HFloat(0.3387183746560611), Theta[5] = HFloat(0.2790527227750356), Theta[6] = HFloat(0.23736308172902912), Theta[7] = HFloat(0.2090446551994883), Theta[8] = HFloat(0.1913243977758092), Theta[9] = HFloat(0.1827101427467962)]

  

[[0, 0], [.5, HFloat(0.3592770908855222)], [1.0, HFloat(0.5657329830840051)], [1.5, HFloat(0.6872174854707988)], [2.0, HFloat(0.76034019145543)], [2.5, HFloat(0.8052695951373906)], [3.0, HFloat(0.8333150468603043)], [3.5, HFloat(0.8508739043275494)], [4.0, HFloat(0.8614824212622565)], [4.5, HFloat(0.8667958074968981)]]

  

[[0, 1], [.5, HFloat(0.7269817315917098)], [1.0, HFloat(0.5464778444851895)], [1.5, HFloat(0.4238106210804488)], [2.0, HFloat(0.3387183746560611)], [2.5, HFloat(0.2790527227750356)], [3.0, HFloat(0.23736308172902912)], [3.5, HFloat(0.2090446551994883)], [4.0, HFloat(0.1913243977758092)], [4.5, HFloat(0.1827101427467962)]]

 (5)
> 

 p3:=plot([data1],color=red,linestyle=1,thickness=3,legend = ["FDM"]):
 p4:=plot([data2],color = red, linestyle = 1, thickness = 3, legend = ["FDM"]):
 display({p1, p3}, axes = boxed);
 display({p2, p4}, axes = boxed);

 

 

NULL

Download FDM2.mw

That dataset...

tomleslie 13836
April 26 2023
0 0

@vs140580 

you are using contains the edge {0, 1}. 0 is not a positive integer, so the type check

Gs::list(set(And(set(posint), 2 &under nops))),

will fail. Suggest you change 'posint' in the above to 'nonnegint', as in the attached.

> 

restart
> 

"GraphsToRGraphs:= proc( Gs::list(set(And(set(nonnegint), 2 &under nops))), file::{string, symbol}:= 'terminal' ) local k, G, E, n, fp, pf:= curry(fprintf, fp); try fp:= `if`(file = 'terminal', file, fopen(file, WRITE)) catch: fp:= 'terminal' end try; pf("[\\n"); for k,G in Gs do E:= [G[]]; n:= nops(E)-1; if k>1 then pf(",\\n\\n") fi; pf(cat(" from <- c(", "%d,"$n, "%d),\\n"), op~(1,E)[]); pf(cat(" to <- c(", "%d,"$n, "%d),\\n"), op~(2,E)[]); pf(" ft <- cbind(from, to),\\n"); pf(" Graph%d <- ftM2graphNEL(ft,edgemode=\"undirected\")", k) od; pf("\\n]\\n"); if file <> 'terminal' then fclose(fp) fi; return end proc : GraphsToRGraphs( [ {{0,1},{1,2},{1,5},{2,3},{2,26},{3,4},{3,19},{4,5},{5,6},{6,7},{7,8},{7,18},{8,9},{9,10},{9,16},{10,11},{10,15},{11,12},{12,13},{13,14},{14,15},{16,17},{17,18},{19,20},{20,21},{20,25},{21,22},{22,23},{23,24},{24,25},{26,27},{27,28},{27,32},{28,29},{29,30},{30,31},{31,32}}, {{1, 2}, {1, 7}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}} ] );"

[
    from <- c(0,1,1,2,2,3,3,4,5,6,7,7,8,9,9,10,10,11,12,13,14,16,17,19,20,20,21,22,23,24,26,27,27,28,29,30,31),
    to <- c(1,2,5,3,26,4,19,5,6,7,8,18,9,10,16,11,15,12,13,14,15,17,18,20,21,25,22,23,24,25,27,28,32,29,30,31,32),
    ft <- cbind(from, to),
    Graph1 <- ftM2graphNEL(ft,edgemode="undirected"),

    from <- c(1,1,2,3,4,5,6),
    to <- c(2,7,3,4,5,6,7),
    ft <- cbind(from, to),
    Graph2 <- ftM2graphNEL(ft,edgemode="undirected")
]

  
> 

NULL

Download notposint.mw

 

If I fix...

tomleslie 13836
April 26 2023
0 0

@KIRAN SAJJAN 

ypur latest worksheet so that it actually executes, I produced the attached. I had to fix loads of syntax and make a lot of guesses aboout your intentions. Somee of these guesses may be wrong.

I also had to use the DirectSearch package from the Maple Application centre in order to solve the reulting algebraic equations. You might be able to get a simple fsolve() to work with these equations by restricting the solution intervals for all variables. I don't have the time/patience to try this

``

> 

restart;
with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, shadebetween, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

 (1)
> 

A := 0.5;
M := 0.5;
R := 0.5;
Pr := 0.72;
Ec := 0.1;
g := 0.1;

.5

  

.5

  

.5

  

.72

  

.1

  

.1

 (2)
> 

ODE := diff(f(eta), eta, eta, eta) + f(eta)*diff(f(eta), eta, eta) - diff(f(eta), eta)^2 - M*diff(f(eta), eta) - A*diff(f(eta), eta) - A*eta*diff(f(eta), eta, eta)/2 = 0, (1 + R)*diff(Theta(eta), eta, eta) + Pr*(f(eta)*diff(Theta(eta), eta) - Theta(eta)*diff(f(eta), eta)) - Pr*A*(Theta(eta) + eta/2*diff(Theta(eta), eta)) + Pr*(Ec*diff(f(eta), eta, eta)^2 + g*Theta(eta));
BCS:=f(0)=0,D(f)(0)=1,Theta(0)=1,D(f)(10)=0,D(Theta)(10)=0;
sol:=dsolve({ODE,BCS},numeric,output = listprocedure):
p1:=odeplot(sol,[eta,f(eta)],0..9,color=green,linestyle=3,thickness=3,legend = ["dsolve"]):
p2:=odeplot(sol, [eta, Theta(eta)], 0 .. 9, color = green, linestyle = 3, thickness = 3, legend = ["dsolve"]):

diff(diff(diff(f(eta), eta), eta), eta)+f(eta)*(diff(diff(f(eta), eta), eta))-(diff(f(eta), eta))^2-1.0*(diff(f(eta), eta))-.2500000000*eta*(diff(diff(f(eta), eta), eta)) = 0, 1.5*(diff(diff(Theta(eta), eta), eta))+.72*f(eta)*(diff(Theta(eta), eta))-.72*Theta(eta)*(diff(f(eta), eta))-.288*Theta(eta)-.1800000000*eta*(diff(Theta(eta), eta))+0.72e-1*(diff(diff(f(eta), eta), eta))^2

  

f(0) = 0, (D(f))(0) = 1, Theta(0) = 1, (D(f))(10) = 0, (D(Theta))(10) = 0

 (3)
> 

d1:= (f[i+1]-f[i-1])/(2*h):
d2:=(f[i+1]-2*f[i]+f[i-1])/h^2:
d3:= (f[i+1]-3*f[i]+3*f[i-1]-f[i-2])/h^3:
d4:= (Theta[i+1]-Theta[i-1])/(2*h):
d5:=(Theta[i+1]-2*Theta[i]+Theta[i-1])/h^2:
h := 0.5:
ode1 := d3 + f[i]*d2 - d1^2 - M*d1 - A*d1 - A*eta*d2/2 = 0:
ode2 := (1 + R)*d5 + Pr*(-d1*Theta[i] + d4*f[i]) - Pr*A*(Theta[i] + eta/2*d4) + Pr*(Ec*d2^2 + g*Theta[i]):
> 

eta[0]:=0:
f[0]:=0:
Theta[0] := 1:
f[-1]:=f[1]-2*h:
f[10]:=f[9]:
Theta[10] := Theta[9]:
for i from 0 to 9 do
    eta[i+1]:=eta[i]+h;
    Eq[i]:=eval(ode1,eta=eta[i]);
    eq[i]:=eval(ode2,eta=eta[i]);
end do:
> 

 soln1:= DirectSearch:-SolveEquations( [Eq[0],Eq[1],Eq[2],Eq[3],Eq[4],Eq[5],Eq[6],Eq[7],Eq[8],Eq[9]])[3];
 soln2:= DirectSearch:-SolveEquations( eval([eq[0],eq[1],eq[2],eq[3],eq[4],eq[5],eq[6],eq[7],eq[8],eq[9]], soln1))[3];
 data1:=[seq([eta[0]+i*h,eval(f[i], soln1 )],i=0..9)]:
 data2 := [seq([h*i + eta[0], eval(Theta[i], soln2)], i = 0 .. 9)]:

[f[-2] = HFloat(-1.8128916235442254), f[1] = HFloat(0.3592770908855222), f[2] = HFloat(0.5657329830840051), f[3] = HFloat(0.6872174854707988), f[4] = HFloat(0.76034019145543), f[5] = HFloat(0.8052695951373906), f[6] = HFloat(0.8333150468603043), f[7] = HFloat(0.8508739043275494), f[8] = HFloat(0.8614824212622565), f[9] = HFloat(0.8667958074968981)]

  

[Theta[-1] = HFloat(1.425809612677366), Theta[1] = HFloat(0.7269817315917098), Theta[2] = HFloat(0.5464778444851895), Theta[3] = HFloat(0.4238106210804488), Theta[4] = HFloat(0.3387183746560611), Theta[5] = HFloat(0.2790527227750356), Theta[6] = HFloat(0.23736308172902912), Theta[7] = HFloat(0.2090446551994883), Theta[8] = HFloat(0.1913243977758092), Theta[9] = HFloat(0.1827101427467962)]

 (4)
> 

 p3:=plot([data1],color=red,linestyle=1,thickness=3,legend = ["FDM"]):
 p4:=plot([data2],color = red, linestyle = 1, thickness = 3, legend = ["FDM"]):
 display({p1, p3}, axes = boxed);
 display({p2, p4}, axes = boxed);

 

 

NULL

Download FDM.mw

Hmmm...

tomleslie 13836
April 24 2023
0 0

@KIRAN SAJJAN 

You state

i want to copare runge kutta 4th order method solution and FEM SOLUTION  comparison 

Basically, you can't! Runge-Kutta methods can only be applied to initial value problems, and you have a boundary value problem. The Maple solution you provide in your latest worksheet, using the command

 
 

msol:=dsolve({subs(beta=0.5,ODE),BCS},numeric,continuation=lambda1):

is actually computed using a finite difference method built-in to Maple. (I have no intention of debugging your attempt at writing your own finite difference method)

And...

tomleslie 13836
April 24 2023
0 0

@KIRAN SAJJAN 

what is the relevance of this latest post? Seems to be nothing to do with your original question?

So precisely what are you trying to achieve?

Don't understand this comment...

tomleslie 13836
April 24 2023
0 0

@KIRAN SAJJAN 

The" shooting mehod" is a way of converting a BVP to an IVP, by appropriate manipulation of the boundary conditions - in general it would not be used just to solve a BVP.

Suggest you read the help at

?dsolve/numeric/BVP

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