tomleslie

If I fix...

Answered: tomleslie 13579

April 26 2023

0 0

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

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

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

Download FDM.mw

Hmmm...

Answered: tomleslie 13579

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

Answered: tomleslie 13579

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

Answered: tomleslie 13579

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

Try...

Answered: tomleslie 13579

April 17 2023

0 0

@lcz

adding the following to my previous answer

 parse~(Explode(Select(IsDigit, sstr1)));

Question...

Answered: tomleslie 13579

April 08 2023

0 0

@vs140580

makes no sense. Try posting a worksheet whicvh illustrates the problem

There are...

Answered: tomleslie 13579

April 06 2023

0 0

@Linhuchong

a variety of ways - but I have to admit that I cheated a little, by using Maple's geometry() package to ectract relevant ellipse parameters, like the centre and lengths of the major and minor axes - see the attached

Download ellgeom.mw

Well...

Answered: tomleslie 13579

April 05 2023

0 2

If you upload your worksheet using the big green up-arrow in the Mapleprimes edit toolbar, then someone can probably help you.

After fixing...

Answered: tomleslie 13579

April 05 2023

0 0

@KIRAN SAJJAN

numerous (minor) typos, the attached produces all the plots you *seem* to want.

PS I threw out a lot of code whichyou don't seem to need.

>

restart: with(plots):

rf := 1050: kf := .52: cpf := 3617: `&sigma:f` := .8:
sigma1 := 4.10*10^7: rs1 := 19300: ks1 := 318: cps1 := 129:
sigma2 := 10^(-10): rs2 := 3970: ks2 := 40: cps2 := 765:
sigma3 := 6.30*10^7: rs3 := 10500: ks3 := 429: cps3 := 235:
sigma4 := 10^(-12): rs4 := 4250: ks4 := 8.9538: cps4 := 686.2:

l := 1: la := 1: Qc := 1: Tp := 1:
g := 1: Ec := 1: H := 3:

Rd := 1: M := 2: Pr := 21: p1 := 0.25e-1: p2 := 0.25e-1: t := (1/3)*Pi: lambda0 := 3: lambda1 := 2: delta1 := 1:
a1 := 1/(1-p1-p2)^2.5:
a1 := 1/(1-p1-p2)^2.5:
a2 := 1-p1-p2+p1*rs1/rf+p2*rs2/rf:
a3 := 1-p1-p2+p1*cps1/cpf+p2*cps2/cpf:
a4 := ((ks1*p1+ks2*p2)/(p1+p2)+2*kf+2*(ks1*p1+ks2*p2)-2*(p1+p2)*kf)/((ks1*p1+ks2*p2)/(p1+p2)+2*kf-ks1*p1-ks2*p2+(p1+p2)*kf):
a5 := 1+3*((p1*sigma1+p2*sigma2)/`&sigma:f`-p1-p2)/(2+(p1*sigma1+p2*sigma2)/((p1+p2)*`&sigma:f`)-((p1*sigma1+p2*sigma2)/`&sigma:f`-p1-p2)):
a6 := 1/(1-p1-p2)^2.5:
a7 := 1-p1-p2+p1*rs3/rf+p2*rs4/rf:
a8 := 1-p1-p2+p1*cps3/cpf+p2*cps4/cpf:
a9 := ((ks3*p1+ks4*p2)/(p1+p2)+2*kf+2*(ks3*p1+ks4*p2)-2*(p1+p2)*kf)/((ks3*p1+ks4*p2)/(p1+p2)+2*kf-ks3*p1-ks4*p2+(p1+p2)*kf):
a10 := 1+3*((p1*sigma3+p2*sigma4)/`&sigma:f`-p1-p2)/(2+(p1*sigma3+p2*sigma4)/((p1+p2)*`&sigma:f`)-((p1*sigma3+p2*sigma4)/`&sigma:f`-p1-p2)):

OdeSys := (a1+6*delta1*(diff(U0(Y), Y))^2)*(diff(U0(Y), Y, Y))-a2*R*(diff(U0(Y), Y))+lambda0-a5*M^2*U0(Y)+la*R*a2*Theta0(Y)*(1+Qc*Theta0(Y)), (a1+g*l+6*delta1*(diff(U0(Y), Y))^2)*(diff(U1(Y), Y, Y))+12*delta1*(diff(U0(Y), Y))*(diff(U1(Y), Y))*(diff(U0(Y), Y, Y))-a2*R*(diff(U1(Y), Y))-(H^2*a2*l+M^2*a5)*U1(Y)+lambda1+la*R*a2*Theta1(Y)*(1+2*Qc*Theta0(Y)), (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y)))*(diff(Theta0(Y), Y, Y))/Pr-a3*R*(diff(Theta0(Y), Y))+a5*Ec*M^2*U0(Y)^2+a1*Ec*(diff(U0(Y), Y))^2+2*Ec*delta1*(diff(U0(Y), Y))^4+4*Rd*(Tp-1)*(diff(Theta0(Y), Y))^2*(1+(2*(Tp-1))*Theta0(Y)+(Tp-1)^2*Theta0(Y)^2)/Pr, (a4+(4*Rd*(1/3))*(1+(Tp-1)^3*Theta0(Y)^3+3*(Tp-1)^2*Theta0(Y)^2+(3*(Tp-1))*Theta0(Y)))*(diff(Theta1(Y), Y, Y))/Pr-a3*R*(diff(Theta1(Y), Y))-l*a3*H^2*Theta1(Y)+2*a5*Ec*M^2*U0(Y)*U1(Y)+2*a1*Ec*(diff(U0(Y), Y))*(diff(U1(Y), Y))+8*Ec*delta1*(diff(U0(Y), Y))^3*(diff(U1(Y), Y))+4*Rd*(Tp-1)*(Theta1(Y)+(2*(Tp-1))*Theta0(Y)*Theta1(Y)+(Tp-1)^2*Theta0(Y)^2*Theta1(Y))*(diff(Theta0(Y), Y, Y))/Pr+8*Rd*(Tp-1)*(diff(Theta0(Y), Y))^2*((Tp-1)^2*Theta0(Y)*Theta1(Y)+(Tp-1)*Theta1(Y))/Pr+8*Rd*(Tp-1)*(diff(Theta0(Y), Y))*(diff(Theta1(Y), Y))*(1+(2*(Tp-1))*Theta0(Y)+(Tp-1)^2*Theta0(Y)^2)/Pr:

Cond := U0(0) = 0, U1(0) = 0, Theta0(0) = 0, Theta1(0) = 0, U0(1) = 0, U1(1) = 0, Theta0(1) = 1, Theta1(1) = 0:
RVals := [1, 2, 3]: Ans1[j] := {}: Ans[j] := {}:

for j to numelems(RVals) do
Ans[j] := dsolve(eval([OdeSys, Cond], R = RVals[j]), numeric, output = listprocedure)
end do:

OdeSys1 := (a6+6*delta1*(diff(V0(Y), Y))^2)*(diff(V0(Y), Y, Y))-a7*R*(diff(V0(Y), Y))+lambda0-a10*M^2*V0(Y)+la*R*a7*Phi0(Y)*(1+Qc*Phi0(Y)), (a6+g*l+6*delta1*(diff(V0(Y), Y))^2)*(diff(V1(Y), Y, Y))+12*delta1*(diff(V0(Y), Y))*(diff(V1(Y), Y))*(diff(V0(Y), Y, Y))-a7*R*(diff(V1(Y), Y))-(H^2*a7*l+M^2*a10)*V1(Y)+lambda1+la*R*a7*Phi1(Y)*(1+2*Qc*Phi0(Y)), (a9+(4*Rd*(1/3))*(1+(Tp-1)^3*Phi0(Y)^3+3*(Tp-1)^2*Phi0(Y)^2+(3*(Tp-1))*Phi0(Y)))*(diff(Phi0(Y), Y, Y))/Pr-a8*R*(diff(Phi0(Y), Y))+a10*Ec*M^2*V0(Y)^2+a6*Ec*(diff(V0(Y), Y))^2+2*Ec*delta1*(diff(V0(Y), Y))^4+4*Rd*(Tp-1)*(diff(Phi0(Y), Y))^2*(1+(2*(Tp-1))*Phi0(Y)+(Tp-1)^2*Phi0(Y)^2)/Pr, (a9+(4*Rd*(1/3))*(1+(Tp-1)^3*Phi0(Y)^3+3*(Tp-1)^2*Phi0(Y)^2+(3*(Tp-1))*Phi0(Y)))*(diff(Phi1(Y), Y, Y))/Pr-a8*R*(diff(Phi1(Y), Y))-l*a8*H^2*Phi1(Y)+2*a10*Ec*M^2*V0(Y)*V1(Y)+2*a6*Ec*(diff(V0(Y), Y))*(diff(V1(Y), Y))+8*Ec*delta1*(diff(V0(Y), Y))^3*(diff(V1(Y), Y))+4*Rd*(Tp-1)*(Phi1(Y)+(2*(Tp-1))*Phi0(Y)*Phi1(Y)+(Tp-1)^2*Phi0(Y)^2*Phi1(Y))*(diff(Phi0(Y), Y, Y))/Pr+8*Rd*(Tp-1)*(diff(Phi0(Y), Y))^2*((Tp-1)^2*Phi0(Y)*Phi1(Y)+(Tp-1)*Phi1(Y))/Pr+8*Rd*(Tp-1)*(diff(Phi0(Y), Y))*(diff(Phi1(Y), Y))*(1+(2*(Tp-1))*Phi0(Y)+(Tp-1)^2*Phi0(Y)^2)/Pr:

Cond1 := V0(0) = 0, V1(0) = 0, Phi0(0) = 0, Phi1(0) = 0, V0(1) = 0, V1(1) = 0, Phi0(1) = 1, Phi1(1) = 0:

for j to numelems(RVals) do
Ans1[j] := dsolve(eval([OdeSys1, Cond1], R = RVals[j]), numeric, output = listprocedure)
end do:

>

  cols := [red, blue, black]:
  plotA:= display
          ( [ seq
              ( odeplot
                ( Ans[k],
                  [Y, U0(Y)],
                  0..1,
                  color=cols[k]
                 ),
                 k=1..numelems(RVals)
               )
             ],
             'axes'= 'boxed',
             labels=[Y,'U0']
            );
  plotB:= display
           ( [ seq
               ( odeplot
                 ( Ans[k],
                   [Y, U1(Y)],
                   0..1,
                   color=cols[k]
                 ),
                 k=1..numelems(RVals)
               )
             ],
             'axes'= 'boxed',
             labels=[Y,'U1']
           );
  plotC:= display
          ( [ seq
              ( odeplot
                ( Ans1[k],
                  [Y, V0(Y)],
                  Y=0..1,
                  color=cols[k]
                ),
                k=1..numelems(RVals)
              )
            ],
           'axes'= 'boxed',
           'linestyle' = 'dashdot',
           labels=[Y,'U0']
         );
  plotD:= display
           ( [ seq
               ( odeplot
                 ( Ans1[k],
                   [Y,V1(Y)],
                   Y=0..1,
                   color=cols[k]
                 ),
                 k=1..numelems(RVals)
               )
             ],
             'axes'= 'boxed',
             'linestyle' = 'dashdot',
             labels=[Y,'U1']
           );
  plotE:= display
          ( [ seq
              ( odeplot
                ( Ans[k],
                  [Y,Theta0(Y)],
                  Y=0..1,
                  color=cols[k]
                ),
                k=1..numelems(RVals)
              )
            ],
            'axes'= 'boxed',
            labels=[Y,'Theta0']
          );
  plotF:= display
           ( [ seq
               ( odeplot
                 ( Ans[k],
                   [Y, Theta1(Y)],
                   Y=0..1,
                   color=cols[k]
                 ),
                 k=1..numelems(RVals)
               )
             ],
             'axes'= 'boxed',
              labels=[Y,'Theta1']
           );
  plotG:= display
          ( [ seq
              ( odeplot
                ( Ans1[k],
                  [Y, Phi0(Y)],
                  Y=0..1,
                  color=cols[k]
                ),
                k=1..numelems(RVals)
              )
           ],
           'axes'= 'boxed',
           'linestyle' = 'dashdot',
           labels=[Y,'Theta0']
         );
  plotH:= display
           ( [ seq
               ( odeplot
                 ( Ans1[k],
                   [Y,Phi1(Y)],
                   Y=0..1,
                   color=cols[k]
                 ),
                 k=1..numelems(RVals)
               )
             ],
             'axes'= 'boxed',
             'linestyle' = 'dashdot',
             labels=[Y,'Theta1']
           );

>	plots:-display([plotE, plotG]); plots:-display([plotF, plotH]);

>

Download odeStuff.mw

Try...

Answered: tomleslie 13579

March 29 2023

0 0

@MANUTTM

op(4, s2)[][];

The function...

Answered: tomleslie 13579

March 24 2023

0 0

'g2' does not contain the name 'q', and so is entirely independent of it.

The expression TC contains the variable names

{Cd, D, O1, O2, P, beta, h1, h2, i, q, r, t1, t2, theta1, theta2}.

However when you differentiate this expression with respect to t2, the result contains only the variables

{Cd, D, O2, beta, h2, i, t2, theta2},

so

{O1, h1, q, t1, theta1}

have been eliminated

fsolve()...

Answered: tomleslie 13579

March 22 2023

1 0

is a numerical solver. Thus every quantity in your equation(s), other than the variable(s) you want fsolve() to compute has to have a numeric value.

I notice that you sayyou are trying to "solve a system of differential equations" for which I would have expected you to be using dsolve(), with, or without the numeric option.

This would be so much easier to diagnose/fix if you uploaded your worksheet using the big green up-arrow in the Mapleprimes toolbar

So why don't you...

Answered: tomleslie 13579

March 20 2023

0 0

just use the built-in help which comes with your Maple release?

As in

Help->What's new

from the menus in the standard GUI

As far as I can tell...

Answered: tomleslie 13579

March 16 2023

0 0

Maple 2017 produces "correct" answer. All subsequent Maple releases produce weird answer.

Using lists rather than expressions sequences appears to make the problem go away in later releases (Maple 2018 - Maple2023)

Can reproduce in Maple2023...

Answered: tomleslie 13579

March 16 2023

0 0

This is a bug, and it is a bad one.

If one uses sequences of names rather than integers, then one can see the difference in the calculations being performed depending on whether the Units() package is loaded or not. The second calculation in the attached makes sense: the first one makes no sense at all.

>

restart;