## 13791 Reputation

14 years, 102 days

## Calculations and some graphs...

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:
 >
 > # First strategy
 > # Second strategy
 > # Third strategy

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)
 > # 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);

## You...

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

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

## Ooooops...

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

 >
 >
 >

using the big green up-arrow in the Mapleprimes toolbar

## A better(?) way...

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];
 (1)
 >

## Before...

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

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

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

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

 > restart; with(plots);
 (1)
 > A := 0.5; M := 0.5; R := 0.5; Pr := 0.72; Ec := 0.1; g := 0.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"]):
 (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)]);
 (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)];
 (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);

## That dataset...

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.

 >
 >
 [     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") ]
 >

## If I fix...

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);
 (1)
 > A := 0.5; M := 0.5; R := 0.5; Pr := 0.72; Ec := 0.1; g := 0.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"]):
 (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)]:
 (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);

## Hmmm...

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

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

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
﻿