Advanced Search

mehdi jafari

769 Reputation

13 Badges

11 years, 334 days
Contact mehdi jafari

MaplePrimes Activity


These are answers submitted by mehdi jafari

check your answer via odetest command...

mehdi jafari 769
April 28 2014
3 1

as maple help page :

The odetest command checks explicit and implicit solutions for ODEs by making a careful simplification of the ODE with respect to the given solution. If the solution is valid, the returned result will be 0; otherwise, the algebraic remaining expression will be returned. In the case of systems of ODEs, odetest can only test explicit solutions, given either as a set or as a list of sets. (For information on non-composed sets of solutions for nonlinear systems, see dsolve,system.)

> 

restart:
> 

sys_ode := {diff(d11(m), m) = -(3*sin(m)^2-1)*d31(m)/a^(3/2)+(-3*cos(m)*sin(m)/a^(3/2))*d41(m), diff(d21(m), m) = (-3*cos(m)*sin(m)/a^(3/2))*d31(m)-(3*cos(m)^2-1)*d41(m)/a^(3/2), diff(d31(m), m) = -a^(3/2)*d11(m), diff(d41(m), m) = -a^(3/2)*d21(m)};

{diff(d11(m), m) = -(3*sin(m)^2-1)*d31(m)/a^(3/2)-3*cos(m)*sin(m)*d41(m)/a^(3/2), diff(d21(m), m) = -3*cos(m)*sin(m)*d31(m)/a^(3/2)-(3*cos(m)^2-1)*d41(m)/a^(3/2), diff(d31(m), m) = -a^(3/2)*d11(m), diff(d41(m), m) = -a^(3/2)*d21(m)}

 (1)
> 

dsolve(sys_ode);

{d11(m) = (1/2)*_C2*sin(2*m)+3*sin(m)*ln(cos(m)+I*sin(m))*_C4-_C3*cos(2*m)-I*cos(m)*_C4+_C1*sin(m), d21(m) = -(1/4)*(-6*ln(cos(m)+I*sin(m))*_C4*sin(4*m)+(6*I)*cos(4*m)*ln(cos(m)+I*sin(m))*_C4+(18*I)*cos(6*m)*ln(cos(m)+I*sin(m))*_C4+(6*I)*cos(2*m)*ln(cos(m)+I*sin(m))*_C4-2*_C1*sin(2*m)-6*_C1*sin(6*m)-2*_C1*sin(4*m)-10*sin(2*m)*_C4+I*_C2*cos(m)-6*ln(cos(m)+I*sin(m))*_C4*sin(2*m)-18*ln(cos(m)+I*sin(m))*_C4*sin(6*m)-3*_C2*sin(7*m)-6*_C4*sin(6*m)+10*_C4*sin(4*m)+5*sin(m)*_C2+2*sin(5*m)*_C2-6*sin(3*m)*_C2+6*_C3*cos(7*m)-2*cos(m)*_C3-4*cos(5*m)*_C3+(6*I)*_C1-(6*I)*_C4+(10*I)*_C3*sin(m)+(3*I)*_C2*cos(7*m)+(6*I)*_C3*sin(7*m)-(4*I)*_C3*sin(5*m)-(2*I)*_C2*cos(5*m)+(2*I)*_C1*cos(4*m)-(10*I)*_C4*cos(4*m)+(6*I)*_C1*cos(6*m)+(6*I)*_C4*cos(6*m)+(6*I)*_C2*cos(3*m)+(2*I)*_C1*cos(2*m)+(10*I)*_C4*cos(2*m)+(18*I)*ln(cos(m)+I*sin(m))*_C4)/(-(3*I)*cos(m)-(3*I)*cos(5*m)+(2*I)*cos(3*m)+3*sin(m)+3*sin(5*m)-2*sin(3*m)), d31(m) = -(1/8)*a^(3/2)*(-12*ln(cos(m)+I*sin(m))*_C4*sin(4*m)-4*_C1*sin(2*m)-12*_C1*sin(6*m)-4*_C1*sin(4*m)+(12*I)*ln(cos(m)+I*sin(m))*_C4*cos(4*m)+(36*I)*ln(cos(m)+I*sin(m))*_C4*cos(6*m)+(12*I)*ln(cos(m)+I*sin(m))*_C4*cos(2*m)+(12*I)*_C1+(24*I)*_C4+(19*I)*_C2*cos(m)+(16*I)*_C2*cos(5*m)+(4*I)*_C1*cos(4*m)+(40*I)*_C4*cos(4*m)+(12*I)*_C1*cos(6*m)-(24*I)*_C4*cos(6*m)-(6*I)*_C2*cos(3*m)+(4*I)*_C1*cos(2*m)-(40*I)*_C4*cos(2*m)+(36*I)*ln(cos(m)+I*sin(m))*_C4+40*sin(2*m)*_C4-12*ln(cos(m)+I*sin(m))*_C4*sin(2*m)-36*ln(cos(m)+I*sin(m))*_C4*sin(6*m)-3*_C2*sin(7*m)+24*_C4*sin(6*m)-40*_C4*sin(4*m)-13*sin(m)*_C2-16*sin(5*m)*_C2+6*sin(3*m)*_C2+6*_C3*cos(7*m)-2*cos(m)*_C3-4*cos(5*m)*_C3+(10*I)*_C3*sin(m)+(3*I)*_C2*cos(7*m)+(6*I)*_C3*sin(7*m)-(4*I)*_C3*sin(5*m))/(-(3*I)*cos(m)-(3*I)*cos(5*m)+(2*I)*cos(3*m)+3*sin(m)+3*sin(5*m)-2*sin(3*m)), d41(m) = (1/8)*a^(3/2)*(-8*_C4*cos(2*m)-8*_C4*cos(4*m)-36*ln(cos(m)+I*sin(m))*_C4-46*_C3*sin(m)-_C2*cos(m)+3*_C2*cos(7*m)+6*_C3*sin(7*m)-40*_C3*sin(5*m)-2*_C2*cos(5*m)-20*_C1*cos(4*m)+12*_C1*cos(6*m)-24*_C4*cos(6*m)+36*_C3*sin(3*m)+20*_C1*cos(2*m)+(36*I)*ln(cos(m)+I*sin(m))*_C4*sin(6*m)-(60*I)*ln(cos(m)+I*sin(m))*_C4*sin(4*m)+(60*I)*ln(cos(m)+I*sin(m))*_C4*sin(2*m)+(20*I)*_C1*sin(2*m)-(8*I)*_C4*sin(2*m)-60*ln(cos(m)+I*sin(m))*_C4*cos(4*m)+36*ln(cos(m)+I*sin(m))*_C4*cos(6*m)+60*ln(cos(m)+I*sin(m))*_C4*cos(2*m)+(5*I)*_C2*sin(m)+(34*I)*_C3*cos(m)+(12*I)*_C1*sin(6*m)-(24*I)*_C4*sin(6*m)-(20*I)*_C1*sin(4*m)-(8*I)*_C4*sin(4*m)-(6*I)*_C3*cos(7*m)+(3*I)*_C2*sin(7*m)-(2*I)*_C2*sin(5*m)+(40*I)*_C3*cos(5*m)-(36*I)*_C3*cos(3*m)-24*_C4-12*_C1)/(-(3*I)*cos(m)-(3*I)*cos(5*m)+(2*I)*cos(3*m)+3*sin(m)+3*sin(5*m)-2*sin(3*m))}

 (2)
> 

odetest(%, sys_ode);

{0}

 (3)
> 

 

 

Download odetest.mw

your problem !...

mehdi jafari 769
April 28 2014
0 0

GenerateMatrix(sys,w);## here is your problem # u can not generate matrix,while some of your elements of w are numbers  !#

what do u exactly want from the output of GenerateMatrix(sys,w) ? maybe we want to check this and solve with another way, please write down what do u expect from the output of that. 

 

> 

restart:
with(LinearAlgebra):
with(plots):
with(PDEtools):
interface(rtablesize=20):

Introduction
 

 

The purpose of this worksheet is to illsutrate how to solve linear elliptic PDEs.  For concreteness, we will focus on the following PDE:

 

"((∂)^(2)u)/(∂x^(2)) +((∂)^(2)u)/(∂y^(2))=f(x,y),"

 

where the source function f(x, y) is given.  If f(x, y) = 0this is known as Laplace's equation, if not it is Poisson's equation.  

The function f is given. We add these boundary condition :

 

u(x,1)=u(x,0)=ux(0,y)=ux(1,y)=0

Stencil
 

In this section, we will develop the solution to the following PDE assuming a stencil for the differential operator involving five points:

> 

pde := diff(u(x,y),x,x) + diff(u(x,y),y,y) - f(x,y);

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y)-f(x, y)

 (2.1)

As usual, we make use of the GenerateStencil procedure:

> 

GenerateStencil := proc(F,N,{orientation:=center,stepsize:=h,showorder:=true,showerror:=false})
    local vars, f, ii, Degri, stencil, Error, unknowns, Indets, ans, Phi, r, n, phi;

 Phi := convert(F, D);
 vars := op(Phi);
n := PDEtools[difforder](Phi);
f := op(1, op(0, Phi));
 if nops([vars]) <> 1 then
r := op(1, op(0, op(0, Phi)))
else r := 1 end if;
phi := f(vars);
 if orientation = center then if type(N, odd) then
ii := [seq(i, i = -(N-1)*(1/2) .. (N-1)*(1/2))]
else ii := [seq(i, i = -N+1 .. N-1, 2)]
end if
elif orientation = left then
ii := [seq(i, i = -N+1 .. 0)]
elif orientation = right then
ii := [seq(i, i = 0 .. N-1)] end if;
stencil := add(a[ii[i]]*subsop(r = op(r, phi)+ii[i]*stepsize, phi), i = 1 .. N);
Error := (D[`$`(r, n)](f))(vars)-stencil;
Error := convert(series(Error, stepsize, N), polynom);
unknowns := seq(a[ii[i]], i = 1 .. N);
Indets := `minus`(`minus`(`minus`(indets(Error), {vars}), {unknowns}), {stepsize});
 Error := collect(Error, Indets, 'distributed');
    ans := solve({coeffs(Error, Indets),unknowns}, parametric = full);
    if (ans=NULL) then:
       print(`Failure: try increasing the number of points in the stencil`);
       return NULL;
    fi:
    stencil := subs(ans,stencil);
    Error := convert(series(`leadterm`(D[r$n](f)(vars) - stencil),stepsize,N+20),polynom);
    Degri := degree(Error,stepsize);
    if (showorder) then:
       print(cat(`This stencil is of order `,Degri));
    fi:
    if (showerror) then:
       print(cat(`This leading order term in the error is `,Error));
    fi:
    convert(D[r$n](f)(vars) = stencil,diff);

end proc:

 

> 

substencil[1] := GenerateStencil(diff(u(x,y),x,x),3);
substencil[2] := GenerateStencil(diff(u(x,y),y,y),3);

`This stencil is of order 0`

  

diff(diff(u(x, y), x), x) = a[-1]*u(x-h, y)+a[0]*u(x, y)+a[1]*u(x+h, y)

  

`This stencil is of order 0`

  

diff(diff(u(x, y), y), y) = a[-1]*u(x, y-h)+a[0]*u(x, y)+a[1]*u(x, y+h)

 (2.2)

Notice that we have assumed the same stepsize h in the x and y directions.  Putting these into PDE gives:

> 

stencil[1] := subs(substencil[1],substencil[2],pde);

a[-1]*u(x-h, y)+2*a[0]*u(x, y)+a[1]*u(x+h, y)+a[-1]*u(x, y-h)+a[1]*u(x, y+h)-f(x, y)

 (2.3)
> 

 

 

aving now obtained a discrete form of the PDE  we now turn our attention to how to exploit it and obtain a numeric solution.  We first need to discretize the domain `in`(x, y, `&x`([0, 1], [0, 1])) over which we seek a solution.  We will assume N+2 lattice points in the x and y directions, respectively.  The coordinates of the lattice point will be explicitly given by x[i] = Z(i) and y[j] = Z(j),

 

> 

Z := i -> (1/(N+1))*i;
x[0] = Z(0),x[N+1] = Z(N+1),y[0] = Z(0),y[N+1] = Z(N+1);

proc (i) options operator, arrow; i/(N+1) end proc

  

x[0] = 0, x[N+1] = 1, y[0] = 0, y[N+1] = 1

 (2.4)
> 

 

Each node is labelled by our approximation to the true solution of the PDE `&approx;`(u[i, j], u(x[i], y[j])).  Boundary conditions will be used to fix u[i, j] at each of the purple boundary nodes, so the goal of the code will be to solve for u[i, j] at each of the white interior nodes.  We also define f[i, j] = f(x[i], y[j]), which allows us to re-write the stencil as:

> 

Subs := seq(seq(u(x+ii*h,y+jj*h)=u[i+ii,j+jj],ii=-1..1),jj=-1..1),f(x,y)=f[i,j]:
stencil[1] := subs(Subs,stencil[1]);
Stencil[1] := unapply(stencil[1],h,i,j,u,f);

a[-1]*u[i, j-1]+a[-1]*u[i-1, j]+2*a[0]*u[i, j]+a[1]*u[i, j+1]+a[1]*u[i+1, j]-f[i, j]

  

proc (h, i, j, u, f) options operator, arrow; a[-1]*u[i-1, j]+a[-1]*u[i, j-1]+2*a[0]*u[i, j]+a[1]*u[i+1, j]+a[1]*u[i, j+1]-f[i, j] end proc

 (2.5)
> 

 

> 

 

> 

## checking your second procedure #
> 

 

> 

N:=1;_f := (x,y) -> 0;

1

  

proc (x, y) options operator, arrow; 0 end proc

 (2.6)
> 

Z := unapply( (1/(N+1))*i,i);
x[0] := Z(0);x[N+1] := Z(N+1);y[0] := Z(0);y[N+1] := Z(N+1);
 h := evalf(Z(1)-Z(0));

proc (i) options operator, arrow; (1/2)*i end proc

  

0

  

1

  

0

  

1

  

.5000000000

 (2.7)
> 

 

> 

for i from 1 to N do:
    # Neumann boundary condition
     u[N+2,i] :=  u[N,i] ;    
     u[-1,i] := u[1,i];
od;
for i from -1 to N+2 do:
     # Dirichlet boundary condition
     u[i,0] := 0;
     u[i,N+1] := 0;
od;

u[1, 1]

  

u[1, 1]

  

0

  

0

  

0

  

0

  

0

  

0

  

0

  

0

  

0

  

0

 (2.8)
> 

f := Array(0..N+1,0..N+1,[seq([seq(evalf(_f(Z(i),Z(j))),i=0..N+1)],j=0..N+1)],datatype=float);

f := Array(0 .. 2, 0 .. 2, {}, datatype = float[8], storage = rectangular, order = Fortran_order)

 (2.9)
> 

f(1..9);

Matrix(1, 9, {(1, 1) = 0., (1, 2) = 0., (1, 3) = 0., (1, 4) = 0., (1, 5) = 0., (1, 6) = 0., (1, 7) = 0., (1, 8) = 0., (1, 9) = 0.})

 (2.10)
> 

sys := [seq(seq(Stencil[1](h,i,j,u,f),i=0..N+1),j=1..N)];

[a[-1]*u[1, 1]+2*a[0]*u[0, 1]+a[1]*u[1, 1], a[-1]*u[0, 1]+2*a[0]*u[1, 1]+a[1]*u[2, 1], a[-1]*u[1, 1]+2*a[0]*u[2, 1]+a[1]*u[1, 1]]

 (2.11)
> 

w := [seq(seq(u[i,j],i=-1..N+2),j=0..N+1)];
  

[0, 0, 0, 0, 0, u[1, 1], u[0, 1], u[1, 1], u[2, 1], u[1, 1], 0, 0, 0, 0, 0]

 (2.12)
> 

GenerateMatrix(sys,w);## here is your problem # u can not generate matrix,while some of your elements of w are zero !#

Error, (in collect) cannot collect 0

  
> 

 

Procedure

 

Using the previous suty in section stencil we can write the procedure to solve the Laplace equation in [0,1]*[0,1] with the boundary condition Neumann conditions on the vertical boundary and Dirichlet boundary condition on the horizontal baoundary. In our study we will use the same stepsize h in x and y direction.

 

> 

PoissonSolve:=proc(N,_f)
local Z,i,h,y,x,sys,w,f,ff,sol,j,u,Data;
# define basic grid parameters
_f:=ff;
Z := unapply( (1/(N+1))*i,i);
x[0] := Z(0);x[N+1] := Z(N+1);y[0] := Z(0);y[N+1] := Z(N+1);
 h := evalf(Z(1)-Z(0));
# Fix the boundary data and the source matrix
for i from 1 to N do:
    # Neumann boundary condition
     u[N+2,i] :=  u[N,i] ;    
     u[-1,i] := u[1,i];
od:
for i from -1 to N+2 do:
     # Dirichlet boundary condition
     u[i,0] := 0;
     u[i,N+1] := 0;
od:
   ff := Array(0..N+1,0..N+1,[seq([seq(evalf(ff(Z(i),Z(j))),i=0..N+1)],j=0..N+1)]);
# Write down the system of equations to solve and solve them
     sys := [seq(seq(Stencil[1](h,i,j,u,ff),i=0..N+1),j=1..N)];
  w := [seq(seq(u[i,j],i=-1..N+2),j=0..N+1)];
  sol := LinearSolve(GenerateMatrix(sys,w));
   # parse the solution vector sol back into "matrix" form
   for i from 0 to N+1 do:
     for j from 0 to N+1 do:
        u[i,j] := sol[j*(N+1)+i]:
     od:
   od:
# generate a 3D plot of the solution using the surfdata command
   Data := [seq([seq([Z(i),Z(j),u[i,j]],i=0..N+1)],j=0..N+1)]:
surfdata(Data,axes=boxed,labels=[`x`,`y`,`u(x,y)`],shading=zhue,style=patchcontour);
end proc:

Here is an example of the output when the source function is set to zero f(x, y) = 0; i.e., when  reduces down to Laplace's equation:

> 

_f := (x,y) -> 0;
PoissonSolve(5,_f);

proc (x, y) options operator, arrow; 0 end proc

  

Error, (in collect) cannot collect 0

  
> 

 

> 

 

 

Download problem2.mw

your code is so vague !...

mehdi jafari 769
April 28 2014
0 2

you always ask a question as complicated as it possible ! and you expect every body to get phd on your code ! 

you should check your code without procedure and by giving value to inputs of your procedure,so that you can find your problem, i just changed some bracket and some problems in your first procedure and your answers changed and i do not know whether yours are correct or mine ! 
thus please check these complicated codes yourself line by line , beacuse you yourself should know which answer is correct and which in not ! and please just ask about programming questions ! good luck !


> 

restart:
with(LinearAlgebra):
with(plots):
with(PDEtools):
interface(rtablesize=20):

Introduction
 

 

The purpose of this worksheet is to illsutrate how to solve linear elliptic PDEs.  For concreteness, we will focus on the following PDE:

 

"((&PartialD;)^(2)u)/(&PartialD;x^(2)) +((&PartialD;)^(2)u)/(&PartialD;y^(2))=f(x,y),"

 

where the source function f(x, y) is given.  If f(x, y) = 0this is known as Laplace's equation, if not it is Poisson's equation.  

The function f is given. We add these boundary condition :

 

u(x,1)=u(x,0)=ux(0,y)=ux(1,y)=0

Stencil
 

In this section, we will develop the solution to the following PDE assuming a stencil for the differential operator involving five points:

> 

pde := diff(u(x,y),x,x) + diff(u(x,y),y,y) - f(x,y);

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y)-f(x, y)

 (2.1)

As usual, we make use of the GenerateStencil procedure:

> 

GenerateStencil := proc(F,N,{orientation:=center,stepsize:=h,showorder:=true,showerror:=false})
    local vars, f, ii, Degri, stencil, Error, unknowns, Indets, ans, Phi, r, n, phi;

 Phi := convert(F, D);
 vars := op(Phi);
n := PDEtools[difforder](Phi);
f := op(1, op(0, Phi));
 if nops([vars]) <> 1 then
r := op(1, op(0, op(0, Phi)))
else r := 1 end if;
phi := f(vars);
 if orientation = center then if type(N, odd) then
ii := [seq(i, i = -(N-1)*(1/2) .. (N-1)*(1/2))]
else ii := [seq(i, i = -N+1 .. N-1, 2)]
end if
elif orientation = left then
ii := [seq(i, i = -N+1 .. 0)]
elif orientation = right then
ii := [seq(i, i = 0 .. N-1)] end if;
stencil := add(a[ii[i]]*subsop(r = op(r, phi)+ii[i]*stepsize, phi), i = 1 .. N);
Error := (D[`$`(r, n)](f))(vars)-stencil;
Error := convert(series(Error, stepsize, N), polynom);
unknowns := seq(a[ii[i]], i = 1 .. N);
Indets := `minus`(`minus`(`minus`(indets(Error), {vars}), {unknowns}), {stepsize});
 Error := collect(Error, Indets, 'distributed');
    ans := solve({coeffs(Error, Indets),unknowns}, parametric = full);
    if (ans=NULL) then:
       print(`Failure: try increasing the number of points in the stencil`);
       return NULL;
    fi:
    stencil := subs(ans,stencil);
    Error := convert(series(`leadterm`(D[r$n](f)(vars) - stencil),stepsize,N+20),polynom);
    Degri := degree(Error,stepsize);
    if (showorder) then:
       print(cat(`This stencil is of order `,Degri));
    fi:
    if (showerror) then:
       print(cat(`This leading order term in the error is `,Error));
    fi:
    convert(D[r$n](f)(vars) = stencil,diff);

end proc:

 

> 

substencil[1] := GenerateStencil(diff(u(x,y),x,x),3);
substencil[2] := GenerateStencil(diff(u(x,y),y,y),3);

`This stencil is of order 0`

  

diff(diff(u(x, y), x), x) = a[-1]*u(x-h, y)+a[0]*u(x, y)+a[1]*u(x+h, y)

  

`This stencil is of order 0`

  

diff(diff(u(x, y), y), y) = a[-1]*u(x, y-h)+a[0]*u(x, y)+a[1]*u(x, y+h)

 (2.2)

Notice that we have assumed the same stepsize h in the x and y directions.  Putting these into PDE gives:

> 

stencil[1] := subs(substencil[1],substencil[2],pde);

a[-1]*u(x-h, y)+2*a[0]*u(x, y)+a[1]*u(x+h, y)+a[-1]*u(x, y-h)+a[1]*u(x, y+h)-f(x, y)

 (2.3)

 

aving now obtained a discrete form of the PDE  we now turn our attention to how to exploit it and obtain a numeric solution.  We first need to discretize the domain `in`(x, y, `&x`([0, 1], [0, 1])) over which we seek a solution.  We will assume N+2 lattice points in the x and y directions, respectively.  The coordinates of the lattice point will be explicitly given by x[i] = Z(i) and y[j] = Z(j),

 

> 

Z := i -> (1/(N+1))*i;
x[0] = Z(0),x[N+1] = Z(N+1),y[0] = Z(0),y[N+1] = Z(N+1);

proc (i) options operator, arrow; i/(N+1) end proc

  

x[0] = 0, x[N+1] = 1, y[0] = 0, y[N+1] = 1

 (2.4)
> 

 

Each node is labelled by our approximation to the true solution of the PDE `&approx;`(u[i, j], u(x[i], y[j])).  Boundary conditions will be used to fix u[i, j] at each of the purple boundary nodes, so the goal of the code will be to solve for u[i, j] at each of the white interior nodes.  We also define f[i, j] = f(x[i], y[j]), which allows us to re-write the stencil as:

> 

Subs := seq(seq(u(x+ii*h,y+jj*h)=u[i+ii,j+jj],ii=-1..1),jj=-1..1),f(x,y)=f[i,j]:
stencil[1] := subs(Subs,stencil[1]);
Stencil[1] := unapply(stencil[1],h,i,j,u,f);

a[-1]*u[i, j-1]+a[-1]*u[i-1, j]+2*a[0]*u[i, j]+a[1]*u[i, j+1]+a[1]*u[i+1, j]-f[i, j]

  

proc (h, i, j, u, f) options operator, arrow; a[-1]*u[i-1, j]+a[-1]*u[i, j-1]+2*a[0]*u[i, j]+a[1]*u[i+1, j]+a[1]*u[i, j+1]-f[i, j] end proc

 (2.5)
> 

 

Procedure

 

Using the previous suty in section stencil we can write the procedure to solve the Laplace equation in [0,1]*[0,1] with the boundary condition Neumann conditions on the vertical boundary and Dirichlet boundary condition on the horizontal baoundary. In our study we will use the same stepsize h in x and y direction.

 

> 

PoissonSolve:=proc(N,_f)
local Z,i,h,y,x,sys,w,f,sol,j,u,Data;
# define basic grid parameters
Z := i -> (1/(N+1))*i;
x[0] = Z(0),x[N+1] = Z(N+1),y[0] = Z(0),y[N+1] = Z(N+1);
 h := evalf(Z(1)-Z(0));
# Fix the boundary data and the source matrix
for i from 1 to N do:
    # Neumann boundary condition
     u[N+2,i] :=  u[N,i] ;    
     u[-1,i] := u[1,i];
od:
for i from -1 to N+2 do:
     # Dirichlet boundary condition
     u[i,0] := 0;
     u[i,N+1] := 0;
od:
   f := Array(0..N+1,0..N+1,[seq([seq(evalf(_f(Z(i),Z(j))),i=0..N+1)],j=0..N+1)],datatype=float);
# Write down the system of equations to solve and solve them
     sys := [seq(seq(Stencil[1](h,i,j,u,f),i=0..N+1),j=1..N)];
  w := [seq(seq(u[i,j],i=-1..N+2),j=0..N+1)];
  sol := LinearSolve(GenerateMatrix(sys,w));
   # parse the solution vector sol back into "matrix" form
   for i from 0 to N+1 do:
     for j from 0 to N+1 do:
        u[i,j] := sol[j*(N+1)+i]:
     od:
   od:
# generate a 3D plot of the solution using the surfdata command
   Data := [seq([seq([Z(i),Z(j),u[i,j]],i=0..N+1)],j=0..N+1)]:
surfdata(Data,axes=boxed,labels=[`x`,`y`,`u(x,y)`],shading=zhue,style=patchcontour);
end proc:

Here is an example of the output when the source function is set to zero f(x, y) = 0; i.e., when  reduces down to Laplace's equation:

> 

f := (x,y) -> 0;
PoissonSolve(10,f);

proc (x, y) options operator, arrow; 0 end proc

  

Error, (in collect) cannot collect 0

  
> 

 

> 

 


Download Question9X.mw

another way,using do loop !...

mehdi jafari 769
April 25 2014
2 0


> 

restart:
> 

a[0]:=1:for n to 10 do a[n]:=a[0]+a[n-1];od:a[4];

5

 (1)
> 

 


Download for_loop.mw

your boundary conditions are incorrect !...

mehdi jafari 769
April 25 2014
0 3

you should change your boundary conditons. u can see here,i changed time range every thing it possibe,but it does not solve and says that your boundary conditions does not match input boundary conditions.


``

> 

restart:
> 

pde:=diff(u(t,x),t$2)=diff(u(t,x),x$2)-sin(u(t,x));
> 

f:=unapply(x^2,x);
> 

IBC := {u(0,x)=f(x),u(t,-50)=0,D[2](u)(t, 50)=0,D[1](u)(0, x)=-diff(f(x),x)};
> 

 

> 

pds := pdsolve(pde, IBC,time=t,range=0..1, numeric,timestep=0.01, spacestep= 0.1);
> 

 

> 

 

diff(diff(u(t, x), t), t) = diff(diff(u(t, x), x), x)-sin(u(t, x))

  

proc (x) options operator, arrow; x^2 end proc

  

{u(0, x) = x^2, u(t, -50) = 0, (D[1](u))(0, x) = -2*x, (D[2](u))(t, 50) = 0}

  

Error, (in pdsolve/numeric/par_hyp) values in range argument must match input boundary conditions

  
> 

p1 := pds:-plot(t = 0);
> 

p2 := pds:-plot(t = 1/10);
> 

p3 := pds:-plot(t = 1/2);
> 

p4 := pds:-plot(t = 1);
> 

p5 := pds:-plot(t = 2);
> 

plots[display]({p1, p2, p3, p4, p5}, title = `Sine Gordaon at t=0,0.1,0.5,1,2`);
> 

 

> 

 

> 

 

``


Download error_with_boundaries.mw

u can also use simplify...

mehdi jafari 769
April 25 2014
0 0

 

NULL

> 

restart
> 

values := ([seq])(x^3+2*x+1, x = 0 .. 10, 0.1e-2):
> 

numelems(%);

10001

 (1)
> 

> 

V := convert(values, Vector)

V := Vector(4, {(1) = ` 1 .. 10001 `*Vector[column], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (2)
> 

M := convert(values, Matrix)

M := Vector(4, {(1) = ` 1 x 10001 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (3)
> 

``

``

 

Download seq.mw

 

solve using assume !...

mehdi jafari 769
April 24 2014
0 1


``

> 

restart
> 

`assuming`([solve(A*cos(x)+I*R*sin(x) = B*exp(-x)+C*exp(x), x)], [x <> 0]);

RootOf(A*cos(_Z)*exp(_Z)+I*R*sin(_Z)*exp(_Z)-B-C*(exp(_Z))^2)

 (1)
> 

`assuming`([solve(A*cos(x)+I*R*sin(x) = B*exp(-x)+C*exp(x))], [x <> 0]);

{A = A, B = A*cos(x)*exp(x)+I*R*sin(x)*exp(x)-C*(exp(x))^2, C = C, R = R, x = x}

 (2)
> 

A := 5;

5

  

6

  

7

 (3)
> 

`assuming`([solve(A*cos(x)+I*R*sin(x) = B*exp(-x)+C*exp(x), x)], [x <> 0, R::integer]);

RootOf(5*cos(_Z)*exp(_Z)+I*R*sin(_Z)*exp(_Z)-6-7*(exp(_Z))^2)

 (4)
> 

A := 5;

5

  

6

  

7

  

9

 (5)
> 

`assuming`([solve(A*cos(x)+I*R*sin(x) = B*exp(-x)+C*exp(x), x)], [x <> 0]);

RootOf(5*cos(_Z)*exp(_Z)+(9*I)*sin(_Z)*exp(_Z)-6-7*(exp(_Z))^2)

 (6)
> 

evalf(%);

.8972759170+2.061869300*I

 (7)
> 

``

``

 

Download equal.mw

your problem...

mehdi jafari 769
April 24 2014
1 0

your problem is that u have set m=4, but your a has three elements, and thus your vector is out of range !

set n=3

 

``

> 

restart
> 

solution_band_matrix := proc (n, a, b, c, d) local X, i, alpha, y, G, z, x, m; m := n; alpha[1] := a[1]; y[1] := c[1]/alpha[1]; z[1] := d[1]/alpha[1]; for i from 2 to m-1 do alpha[i] := a[i]-b[i]*y[i-1]; y[i] := c[i]/alpha[i]; z[i] := (d[i]-b[i]*z[i-1])/alpha[i] end do; alpha[m] := a[m]-b[m]*y[m-1]; z[m] := (d[m]-b[m]*z[m-1])/alpha[m]; x[m] := z[m]; for i to m-1 do x[i] := z[i]-y[i]*x[i+1] end do end proc:
> 

a := [3, 3, 3];

[3, 3, 3]

  

[0, 1, 1]

  

[1, 1, 0]

  

[0, 1, 1]

  

3

 (1)
> 

solution_band_matrix(n, a, b, c, d);

2/7

 (2)
> 

``

``

Download problem.mw

one way, using Equate...

mehdi jafari 769
April 24 2014
1 2

u can do sth like this :

example.mw

restart:
A:=Matrix(3,[[1,2,3],[4,8,6],[7,8,9]]);
X:=Vector(3,symbol=x);
F:=Vector(3,[8,8,9]);
Equate(X,A^(-1).F);

after closing the program,everything res...

mehdi jafari 769
April 23 2014
0 0

Each time you close maple software,it clears internal memory,just like you use restart command in maple.
as maple help page :
After restart, the variable assignment, loaded package, interface verboseproc variable, and typesetting rules are reset.  The interface elision and typesetting variables remain.for more information,see ?restart

for example :

 

> 

restart:
> 

f:=52;

52

 (1)
> 

f;

52

 (2)
> 

interface(typesetting = standard);

standard

 (3)
> 

restart:
> 

f;

f

 (4)
> 

interface(typesetting);

standard

 (5)
> 

 

 

Download restart.mw

what do u want exactly?...

mehdi jafari 769
April 22 2014
1 1

what do u want exactly? u can use select if  u want to select expressions including something.maybe u want to see these links also :

http://www.mapleprimes.com/questions/148812-How-To-Factor-X-In-This-Expression-

http://www.mapleprimes.com/questions/151261-Problem-With-Collect


> 

restart:
> 

eq:=(1/20)*t[1]^2*(diff(phi(x), x, x))^2*((t[1]+t[2])^5-t[1]^5)/(t[2]^2*E[T])+(1/8)*t[1]^2*(diff(phi(x), x, x))^2*(-2*t[1]-2*t[2])*((t[1]+t[2])^4-t[1]^4)/(t[2]^2*E[T])+(1/3*((1/4)*t[1]^2*(diff(phi(x), x, x))^2*(t[1]+t[2])^2/(t[2]^2*E[T])+(1/4)*t[1]^2*(diff(phi(x), x, x))^2*(-2*t[1]-2*t[2])^2/(t[2]^2*E[T])+(1/2)*t[1]*(diff(phi(x), x, x))*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])/t[2]-t[1]*(diff(psi(x), x))*(-(cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))/t[2]-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))/t[2])/t[2]+(1/2)*t[1]^2*eta(x)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))/t[2]^2+(1/2)*t[1]^2*phi(x)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))/t[2]^2-t[1]*(diff(phi(x), x))*(-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))/t[2]-(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))/t[2])/t[2]+t[1]^2*psi(x)*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*(diff(phi(x), x, x))/t[2]^2))*((t[1]+t[2])^3-t[1]^3)+(1/2*(t[1]^2*psi(x)*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])/t[2]^2+t[1]*(diff(phi(x), x))*(t[1]+t[2])*(-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))/t[2]-(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))/t[2])/t[2]-t[1]*(diff(phi(x), x))*(cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])/t[2]+t[1]*(diff(psi(x), x))*(t[1]+t[2])*(-(cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))/t[2]-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))/t[2])/t[2]-t[1]*(diff(psi(x), x))*((cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])/t[2]+(1/2)*t[1]^2*eta(x)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])/t[2]^2+(1/2)*t[1]^2*phi(x)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])/t[2]^2+(1/4)*t[1]^2*(diff(phi(x), x, x))^2*(t[1]+t[2])^2*(-2*t[1]-2*t[2])/(t[2]^2*E[T])+(1/2)*t[1]*(diff(phi(x), x, x))*(-2*t[1]-2*t[2])*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])/t[2]))*((t[1]+t[2])^2-t[1]^2)+t[1]*phi(x)*((cos(a)^4/E[A]+cos(a)^2*sin(a)^2*(-2*nu[A]/E[A]+1/G[A])+sin(a)^4/E[T])*t[1]*phi(x)/t[2]+(sin(a)^2*cos(a)^2*(1/E[A]+1/E[T]-1/G[A])-(sin(a)^4+cos(a)^4)*nu[A]/E[A])*t[1]*eta(x)/t[2]+(1/2)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/t[2]+sin(a)*cos(a)*(cos(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+sin(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*psi(x)/t[2])+t[1]*(diff(phi(x), x))*(t[1]+t[2])*(cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])+t[1]*eta(x)*((sin(a)^2*cos(a)^2*(1/E[A]+1/E[T]-1/G[A])-(sin(a)^4+cos(a)^4)*nu[A]/E[A])*t[1]*phi(x)/t[2]+(sin(a)^4/E[A]+cos(a)^2*sin(a)^2*(-2*nu[A]/E[A]+1/G[A])+cos(a)^4/E[T])*t[1]*eta(x)/t[2]+(1/2)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/t[2]+sin(a)*cos(a)*(sin(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+cos(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*psi(x)/t[2])+t[1]*psi(x)*(sin(a)*cos(a)*(cos(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+sin(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*phi(x)/t[2]+sin(a)*cos(a)*(sin(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+cos(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*eta(x)/t[2]+cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*(diff(phi(x), x, x))*(t[1]+t[2])^2/t[2]+(4*sin(a)^2*cos(a)^2*(1/E[A]+2*nu[A]/E[A]+1/E[T])+(sin(a)^2-cos(a)^2)^2/G[A])*t[1]*psi(x)/t[2])+t[1]*(diff(psi(x), x))*(t[1]+t[2])*((cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2]):
> 

select(has,eq,sin(a));

(1/3)*((1/4)*t[1]^2*(diff(diff(phi(x), x), x))^2*(t[1]+t[2])^2/(t[2]^2*E[T])+(1/4)*t[1]^2*(diff(diff(phi(x), x), x))^2*(-2*t[1]-2*t[2])^2/(t[2]^2*E[T])+(1/2)*t[1]*(diff(diff(phi(x), x), x))*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])/t[2]-t[1]*(diff(psi(x), x))*(-(cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))/t[2]-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))/t[2])/t[2]+(1/2)*t[1]^2*eta(x)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*(diff(diff(phi(x), x), x))/t[2]^2+(1/2)*t[1]^2*phi(x)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*(diff(diff(phi(x), x), x))/t[2]^2-t[1]*(diff(phi(x), x))*(-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))/t[2]-(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))/t[2])/t[2]+t[1]^2*psi(x)*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*(diff(diff(phi(x), x), x))/t[2]^2)*((t[1]+t[2])^3-t[1]^3)+(1/2)*(t[1]^2*psi(x)*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*(diff(diff(phi(x), x), x))*(-2*t[1]-2*t[2])/t[2]^2+t[1]*(diff(phi(x), x))*(t[1]+t[2])*(-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))/t[2]-(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))/t[2])/t[2]-t[1]*(diff(phi(x), x))*(cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])/t[2]+t[1]*(diff(psi(x), x))*(t[1]+t[2])*(-(cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))/t[2]-cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))/t[2])/t[2]-t[1]*(diff(psi(x), x))*((cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])/t[2]+(1/2)*t[1]^2*eta(x)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*(diff(diff(phi(x), x), x))*(-2*t[1]-2*t[2])/t[2]^2+(1/2)*t[1]^2*phi(x)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*(diff(diff(phi(x), x), x))*(-2*t[1]-2*t[2])/t[2]^2+(1/4)*t[1]^2*(diff(diff(phi(x), x), x))^2*(t[1]+t[2])^2*(-2*t[1]-2*t[2])/(t[2]^2*E[T])+(1/2)*t[1]*(diff(diff(phi(x), x), x))*(-2*t[1]-2*t[2])*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])/t[2])*((t[1]+t[2])^2-t[1]^2)+t[1]*phi(x)*((cos(a)^4/E[A]+cos(a)^2*sin(a)^2*(-2*nu[A]/E[A]+1/G[A])+sin(a)^4/E[T])*t[1]*phi(x)/t[2]+(sin(a)^2*cos(a)^2*(1/E[A]+1/E[T]-1/G[A])-(sin(a)^4+cos(a)^4)*nu[A]/E[A])*t[1]*eta(x)/t[2]+(1/2)*(-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2/t[2]+sin(a)*cos(a)*(cos(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+sin(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*psi(x)/t[2])+t[1]*(diff(phi(x), x))*(t[1]+t[2])*(cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+(sin(a)^2/G[T]+cos(a)^2/G[A])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])+(1/2)*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2*((-cos(a)^2*nu[A]/E[A]-sin(a)^2*nu[T]/E[T])*t[1]*phi(x)/t[2]+(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*eta(x)/t[2]+(1/2)*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2/(E[T]*t[2])+2*cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*psi(x)/t[2])+t[1]*eta(x)*((sin(a)^2*cos(a)^2*(1/E[A]+1/E[T]-1/G[A])-(sin(a)^4+cos(a)^4)*nu[A]/E[A])*t[1]*phi(x)/t[2]+(sin(a)^4/E[A]+cos(a)^2*sin(a)^2*(-2*nu[A]/E[A]+1/G[A])+cos(a)^4/E[T])*t[1]*eta(x)/t[2]+(1/2)*(-sin(a)^2*nu[A]/E[A]-cos(a)^2*nu[T]/E[T])*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2/t[2]+sin(a)*cos(a)*(sin(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+cos(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*psi(x)/t[2])+t[1]*psi(x)*(sin(a)*cos(a)*(cos(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+sin(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*phi(x)/t[2]+sin(a)*cos(a)*(sin(a)^2*(2/E[A]+2*nu[A]/E[A]-1/G[A])+cos(a)^2*(-2*nu[A]/E[A]-2/E[T]+1/G[A]))*t[1]*eta(x)/t[2]+cos(a)*sin(a)*(-nu[A]/E[A]+nu[T]/E[T])*t[1]*(diff(diff(phi(x), x), x))*(t[1]+t[2])^2/t[2]+(4*sin(a)^2*cos(a)^2*(1/E[A]+2*nu[A]/E[A]+1/E[T])+(sin(a)^2-cos(a)^2)^2/G[A])*t[1]*psi(x)/t[2])+t[1]*(diff(psi(x), x))*(t[1]+t[2])*((cos(a)^2/G[T]+sin(a)^2/G[A])*t[1]*(diff(psi(x), x))*(t[1]+t[2])/t[2]+cos(a)*sin(a)*(1/G[A]-1/G[T])*t[1]*(diff(phi(x), x))*(t[1]+t[2])/t[2])

 (1)
> 

 

 

Download select.mw

import your data to a .txt file,and read...

mehdi jafari 769
April 20 2014
0 0

import your data to a .txt file or save them in a .txt file,and read them like this :

restart:
A:=Vector(15):
FileTools[Text][Open]( "testfile.txt" );
0
for i to 15 do
A[i]:=FileTools[Text][ReadFloat]( "testfile.txt" );
od:

if testfile.txt contains 15 number of data .

good luck



set an initial point...

mehdi jafari 769
April 20 2014
2 0

restart:with(Optimization):

Minimize(x^2 + y^2 + 25*(sin(x)^2+sin(y)^2), x=-2*Pi .. 2*Pi , y= -2*Pi .. 2*Pi,initialpoint = {x =Pi}, assume = nonnegative);

f:= (x,y)->cos(x)*sin(y) -(x/(y^2+1));

Maximize((f(x,y), x = -1 .. 2, y = -1 .. 1));

Minimize((f(x,y), x = -1 .. 2, y = -1 .. 1));

as maple help page :

The Minimize and Maximize commands use various methods implemented in a built-in library provided by the Numerical Algorithms Group (NAG). The solvers are iterative in nature and require an initial point. The quality of the solution can depend greatly on the point chosen, especially for nonlinear problems. It is recommended that you provide a point using the initialpoint option. Otherwise, a point is automatically generated.

good luck!

some problems...

mehdi jafari 769
April 18 2014
0 1

1.an unbalanced parantheses in first limit !

2. make some paramteres local

3. your guess1 depends on itself ! thus it is in infinite loop and it can not be calculated ! change your guess 1 !

restart:
fp:=diff(x^2+1,x);

fp2:=diff(x^2+1,x)-1;guess0:=0.1;

guess1:=evalf(guess0-subs(x=guess0,fp)/subs(x=guess0,limit((fp-fp2)/guess0-guess1,x=guess1)));
## it returns this error : 

Error, recursive assignment # 

this is your corrected 1 and 2 problem, but please correct the thrid !

restart:
Secant:=proc(f,x0,tol) local guess1,count,guess0;
fp=diff(f,x);
fp2=diff(f,x)-1;
guess0=x0;
guess1:=evalf(guess0-subs(x=guess0,fp)/subs(x=guess0,limit((fp-fp2)/guess0-guess1,x=guess1))); count:=1;
while abs(guess0-guess1)>tol do
guess0:=guess1;
guess1:=evalf(guess0-subs(x=guess0,fp)/subs(x=guess0,limit((fp-fp2)/guess0-guess1),x=guess1)); count:=count+1;
end;
guess1;count;
end:

maybe it is maybe it is not ! or another...

mehdi jafari 769
April 18 2014
0 1

restart:

int(exp(-z^2*sin(z)^2), z = 0 .. infinity, numeric, epsilon = 0.1e-1); #which returns 2.835068335#

exp(-z^2*sin(z)^2);

F:=convert(%,trig);convert(F,exp);simplify(convert(%,trig));

int(convert(F,exp), z = 0 .. infinity,numeric, epsilon = 0.1e-1); # which returns 2.259350091#

i do not know why the fire can not be uploaded :| . what's up in mapleprimes ?

 

3 4 5 6 7 8 9 Page 5 of 9