tomleslie

13791 Reputation

19 Badges

14 years, 102 days

MaplePrimes Activity


These are replies submitted by tomleslie

@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

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'

@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

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

using the big green up-arrow in the Mapleprimes toolbar

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

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

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

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

 

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

@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

 

@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

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

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

@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