Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@vv I had not expected the maxsols option to work in this case since according to the documentation, "this option is valid only for a single univariate polynomial equation" [my emphasis].  Perhaps the documentation needs to be updated.

@Neel Your hand-written notes on solving the forced wave equation are on the right track.  I have written my take on your notes and have completed the calculations in the following worksheet.  Perhaps this can be of some help.

I must point out that very little of Maple's computational power is used in this worksheet.  These calculations could have been easily done by hand.  Don't get distracted by some of the obscure Maple commands that you may encounter in the worksheet..  Instead, focus on what they are meant to do.

Series solution of the periodically forced wave equation

Here I go step-by-step through the process of obtaining a series solution of the periodically

forced wave equation.

 

Rouben Rostamian

2020-07-08

restart;

Here is the periodically forced wave equation:

pde := diff(w(x,t),t,t) = c^2*diff(w(x,t),x,x) + kappa(x)*cos(Omega*t);

diff(diff(w(x, t), t), t) = c^2*(diff(diff(w(x, t), x), x))+kappa(x)*cos(Omega*t)

We solve it subject to the boundary conditions u(0, t) = 0, u(L, t) = 0 and initial conditions u(x, 0) = phi(x)
and u__t(x, 0) = psi(x).

Following the traditional approach, first we solve the force-free problem and obtain the so-called

homogeneous solution w__h, and then find a particular solution w__p to the forced problem.  Then the

solution will be w = w__h+w__p.

 

The homogeneous problem

 

The homogeneous PDE is

pde_h := diff(w(x,t),t,t) = c^2*diff(w(x,t),x,x);

diff(diff(w(x, t), t), t) = c^2*(diff(diff(w(x, t), x), x))

We look for separable solutions of the form w(x, t) = X(x)*T(t).  Plugging that

into the PDE we get

eval(pde_h, w(x,t)=X(x)*T(t));
pde_separated := % / (X(x)*T(t));   # separate the variables

X(x)*(diff(diff(T(t), t), t)) = c^2*(diff(diff(X(x), x), x))*T(t)

(diff(diff(T(t), t), t))/T(t) = c^2*(diff(diff(X(x), x), x))/X(x)

One side of that equation is a function of t, the other side is a function of x.

But x and t are independent variables, therefore each side has better be a constant.

Thus, we arrive at two differential equations

eq_T := lhs(pde_separated) = -omega^2;
eq_X := rhs(pde_separated) = -omega^2;

(diff(diff(T(t), t), t))/T(t) = -omega^2

c^2*(diff(diff(X(x), x), x))/X(x) = -omega^2

We solve the two ODEs:

sol_T := dsolve(eq_T);

T(t) = _C1*sin(omega*t)+_C2*cos(omega*t)

sol_X_temporary1 := dsolve(eq_X);

X(x) = _C1*sin(omega*x/c)+_C2*cos(omega*x/c)

The expression omega*x/c that occurs above is not very convenient.  We would much rather

have lambda*x instead.  For that, we define lambda in terms of omega through the equation omega*x/c = lambda*x,

whence lambda = omega/c, or equivalently, omega = lambda*c.  In terms of lambda,  sol_X_temporary1 takes the form

sol_X_temporary2 := eval(sol_X_temporary1, omega=lambda*c);

X(x) = _C1*sin(lambda*x)+_C2*cos(lambda*x)

The X solution should satisfy the boundary conditions, that is, X(0) = 0.  Plugging x = 0 in

the solution, we get

eval(sol_X_temporary2, x=0);

X(0) = _C2

and therefore _C2 = 0 and the solution reduces to _C1*sin(lambda*x).  The _C1 factor is

immaterial since any constant multiple of a solution of pde__h is also a solution.
Thus, we see that X(x) is given by

sol_X := subs(_C1=1, _C2=0, sol_X_temporary2);

X(x) = sin(lambda*x)

The solution should also satisfy the boundary conditions X(L) = 0.  Plugging x = L in

the expression for solution we get

eval(sol_X, x=L);

X(L) = sin(lambda*L)

whence lambda*L = n*Pi for any integer "n."  Let's write the resulting lambda as `λ__n` = n*Pi/L, and, recalling

the definition of "lambda,"also set `ω__n` = `λ__n`*c and `λ__n`*c = n*Pi*c/L.

We conclude that there are infinitely many solution for "X ."`` Call them X__n(x) = sin(`λ__n`*x), n = 1, 2, () .. ()

We don't care for non-positive values of n since replacing n by -n is equivalent to multiplying X__n by -1.

 

The functions X__n form a complete orthogonal set under the usual inner product on the

space of square-integrable functions on the interval 0, L.   Let's verify the orthogonality:

Int(sin(lambda[i]*x)*sin(lambda[j]*x), x=0..L);
eval(%, {lambda[i]=i*Pi/L, lambda[j]=j*Pi/L}):
value(%) assuming i::integer, j::integer;

Int(sin(lambda[i]*x)*sin(lambda[j]*x), x = 0 .. L)

0

The functions don't have unit norm.  We could insert a factor of sqrt(2/L) in the definition

of X__n to make them into an orthonormal family, but we won't bother.

Int(sin(lambda[i]*x)^2, x=0..L);
eval(%, lambda[i]=i*Pi/L):
value(%) assuming i::integer;

Int(sin(lambda[i]*x)^2, x = 0 .. L)

(1/2)*L

In view of the solution for T obtained earlier in sol_T, we express the general solution of the homogeneous problem as an infinite sum, where the coefficients a__n and b__n are arbitrary constants.

w__h := Sum((a[n]*cos(omega[n]*t) + b[n]*sin(omega[n]*t))*X[n](x), n=1..infinity);

Sum((a[n]*cos(omega[n]*t)+b[n]*sin(omega[n]*t))*X[n](x), n = 1 .. infinity)

 

 

Orthogonal expansion

 

As has already been noted, the family of functions X__n, n = 1, 2, () .. () form a complete orthogonal set on the space

square-integrable functions.  Any such function, say "f(x)," may be written as an infinite linear combination

of the X__ns, as in
"f(x)=(∑)`alpha__n`*`X__n`(x)."

The coefficients `α__n` are easily determined due to the orthogonality property.  To see that,

multiply the expression above by X__m(x) for an unspecified index m, and integrate the

result over the interval
0, L; int(f(x)*X__m(x), x = 0 .. L) = int(sum(`α__n`*X__n(x)*X__m(x), n = 1 .. infinity), x = 0 .. L) and int(sum(`α__n`*X__n(x)*X__m(x), n = 1 .. infinity), x = 0 .. L) = sum(`α__n`*(int(X__n(x)*X__m(x), x = 0 .. L)), n = 1 .. infinity)
.

Each of the infinitely many integrals on the right-hand side is zero, except for the one
corresponding to n = m which equals "L/(2)."  Therefore the expression above reduces to
"(∫)[0]^(L)f(x)*`X__m`(x) ⅆx=(L)/(2)*`alpha__m` ,"
whence
`α__m` = 2*(int(f(x)*X__m(x), x = 0 .. L))/L

 

The particular solution

 

We have taken the forcing term to be kappa(x)*cos(Omega*t).   We expand the coefficient kappa(x)

in terms of the eigenfunctions X__n  as described in the previous section.  We get
kappa(x) = sum(k__n*X__n(x), n = 1 .. infinity),  where  "`k__n`=2/(L)*(∫)[0]^(L)kappa(x)*`X__n`(x) ⅆx." 

While we are at it, we also expand the initial conditions phi(x) and psi(x) in the same way:

phi(x) = sum(A__n*X__n(x), n = 1 .. infinity),  where  
"`A__n`=2/(L)*(∫)[0]^(L)phi(x)*`X__n`(x) ⅆx,  "
psi(x) = sum(B__n*X__n(x), n = 1 .. infinity),  where  
"`B__n`=2/(L)*(∫)[0]^(L)psi(x)*`X__n`(x) ⅆx.  "

We look for a particular solution of the form W(x)*cos(Omega*t).  Express the yet unknown function
W(x) also as an infinite sum of the eigenfunctions, that is
W(x) = sum(h__n*X__n(x), n = 1 .. infinity),
where h__n are to be determined.

 

Now let us compute.  Form of the sought after solution is:

w_form := w(x,t) = Sum(h[n]*X[n](x), n=1..infinity)*cos(Omega*t);

w(x, t) = (Sum(h[n]*X[n](x), n = 1 .. infinity))*cos(Omega*t)

and the form of the forcing term is (Recall that that the coefficients k__n are known/computable)

force_form := kappa(x) = Sum(k[n]*X[n](x), n=1..infinity);

kappa(x) = Sum(k[n]*X[n](x), n = 1 .. infinity)

eval(pde, [w_form, force_form]);    # plug the above into the PDE
lhs(%) - rhs(%) = 0;                # bring all terms to the left
tmp1 := simplify(%/cos(Omega*t));    # divide through by the common factor cos(Omega*t)

-(Sum(h[n]*X[n](x), n = 1 .. infinity))*Omega^2*cos(Omega*t) = c^2*(Sum(h[n]*(diff(diff(X[n](x), x), x)), n = 1 .. infinity))*cos(Omega*t)+(Sum(k[n]*X[n](x), n = 1 .. infinity))*cos(Omega*t)

-(Sum(h[n]*X[n](x), n = 1 .. infinity))*Omega^2*cos(Omega*t)-c^2*(Sum(h[n]*(diff(diff(X[n](x), x), x)), n = 1 .. infinity))*cos(Omega*t)-(Sum(k[n]*X[n](x), n = 1 .. infinity))*cos(Omega*t) = 0

-(Sum(h[n]*X[n](x), n = 1 .. infinity))*Omega^2-c^2*(Sum(h[n]*(diff(diff(X[n](x), x), x)), n = 1 .. infinity))-(Sum(k[n]*X[n](x), n = 1 .. infinity)) = 0

Recall, however, that X__n is defined as the solution of the differential equation diff(X__n(x), x, x) = -`ω__n`^2*X__n(x)/c^2, therefore

we substitute for it in the previous equation, and simplify:

eval(tmp1, diff(X[n](x),x,x) = -omega[n]^2/c^2*X[n](x) );
combine(%);NULL
tmp2 := factor(%);

-(Sum(h[n]*X[n](x), n = 1 .. infinity))*Omega^2-c^2*(Sum(-h[n]*omega[n]^2*X[n](x)/c^2, n = 1 .. infinity))-(Sum(k[n]*X[n](x), n = 1 .. infinity)) = 0

Sum(h[n]*X[n](x)*omega[n]^2-Omega^2*h[n]*X[n](x)-k[n]*X[n](x), n = 1 .. infinity) = 0

Sum(-X[n](x)*(Omega^2*h[n]-h[n]*omega[n]^2+k[n]), n = 1 .. infinity) = 0

In that last equation we have an infinite linear combination of the basis functions X__n , equaled to zero.

It follows (from orthogonality) that each of the coefficients is zero.  So we set the coefficients equal

to zero and from there compute `h__n:`

op(1, lhs(tmp2)):
coeff(%, X[n](x)) = 0;
h_form := isolate(%, h[n]);

-Omega^2*h[n]+h[n]*omega[n]^2-k[n] = 0

h[n] = k[n]/(-Omega^2+omega[n]^2)

Finally, we substitute the value of h__n computed above in the w_form to obtain the

sought-after particular solution:

w__p := subs(h_form, rhs(w_form));

(Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2), n = 1 .. infinity))*cos(Omega*t)

This is good provided that the forcing frequency Omega avoids the natural frequencies `ω__n`.

 

The solution

 

The solution of the original initial value is the sum of the homogeneous and particular solutions

calculated in the previous sections:

w_sol_general := w__h + w__p;

Sum((a[n]*cos(omega[n]*t)+b[n]*sin(omega[n]*t))*X[n](x), n = 1 .. infinity)+(Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2), n = 1 .. infinity))*cos(Omega*t)

Our next task is to determine the coefficients a__n and b__n so that the initial conditions are satisfied.

Here we will appeal to the infinite series expansions of the initial conditions phi(x) and psi(x) obtained

earlier.  We want w(x, 0) = phi(x), therefore

eval(w_sol_general, t=0) = Sum(A[n]*X[n](x), n=1..infinity);
lhs(%) - rhs(%) = 0;
tmp2 := combine(%);

Sum(a[n]*X[n](x), n = 1 .. infinity)+Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2), n = 1 .. infinity) = Sum(A[n]*X[n](x), n = 1 .. infinity)

Sum(a[n]*X[n](x), n = 1 .. infinity)+Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2), n = 1 .. infinity)-(Sum(A[n]*X[n](x), n = 1 .. infinity)) = 0

Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2)-A[n]*X[n](x)+a[n]*X[n](x), n = 1 .. infinity) = 0

From the orthogonality of the X__n it follows that each of the terms in the sum is zero,

whence we calculate a__n:

op(1, lhs(tmp2));
isolate(%, a[n]);
a_form := expand(%);

k[n]*X[n](x)/(-Omega^2+omega[n]^2)-A[n]*X[n](x)+a[n]*X[n](x)

a[n] = (-k[n]*X[n](x)/(-Omega^2+omega[n]^2)+A[n]*X[n](x))/X[n](x)

a[n] = -k[n]/(-Omega^2+omega[n]^2)+A[n]

Furthermore, we want w__t(x, 0) = psi(x), therefore

diff(w_sol_general, t);
eval(%, t=0) = Sum(B[n]*X[n](x), n=1..infinity);
lhs(%) - rhs(%) = 0;
tmp3 := combine(%);

Sum((-a[n]*omega[n]*sin(omega[n]*t)+b[n]*omega[n]*cos(omega[n]*t))*X[n](x), n = 1 .. infinity)-(Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2), n = 1 .. infinity))*Omega*sin(Omega*t)

Sum(b[n]*omega[n]*X[n](x), n = 1 .. infinity) = Sum(B[n]*X[n](x), n = 1 .. infinity)

Sum(b[n]*omega[n]*X[n](x), n = 1 .. infinity)-(Sum(B[n]*X[n](x), n = 1 .. infinity)) = 0

Sum(b[n]*omega[n]*X[n](x)-B[n]*X[n](x), n = 1 .. infinity) = 0

It follows that

op(1, lhs(tmp3));
b_form := isolate(%, b[n]);

b[n]*omega[n]*X[n](x)-B[n]*X[n](x)

b[n] = B[n]/omega[n]

Finally, we substitute the expressions for a__n and b__n  in w_sol_general and arrive

at the solution of the problem:

w(x,t) = subs(a_form, b_form, w_sol_general);

w(x, t) = Sum(((-k[n]/(-Omega^2+omega[n]^2)+A[n])*cos(omega[n]*t)+B[n]*sin(omega[n]*t)/omega[n])*X[n](x), n = 1 .. infinity)+(Sum(k[n]*X[n](x)/(-Omega^2+omega[n]^2), n = 1 .. infinity))*cos(Omega*t)

The two summations above may be combined into one (although that's not strictly necessary).  It's difficult

doing that in Maple, so here I will do it by hand:

w(x,t) = Sum(((-k[n]/(-Omega^2 + omega[n]^2) + A[n])*cos(omega[n]*t) + B[n]*sin(omega[n]*t)/omega[n] + k[n]*cos(Omega*t)/(-Omega^2 + omega[n]^2))*X[n](x), n=1..infinity);

w(x, t) = Sum(((-k[n]/(-Omega^2+omega[n]^2)+A[n])*cos(omega[n]*t)+B[n]*sin(omega[n]*t)/omega[n]+k[n]*cos(Omega*t)/(-Omega^2+omega[n]^2))*X[n](x), n = 1 .. infinity)

 

 


 

Download wave-equation-forced-vibration.mw

 

After deriving equation (22), the authors identify three special cases for which solutions are available in the literature.  Let's look at their Case 1, that is, a/L=0 and mass_ratio=0.  The equation then should reduce, and they assert that it does, to the classical clamped-free cantilever model.  But I don't see that.

Setting a/L=0 and mass_ratio=0 in equation (22), I get 0=0.  (I did that by hand; no need for Maple for that.)  But that's not correct.  The equation should reduce to cos(βL) cosh(βL) = −1.

Can you verify my calculation?  If I am correct, then there must be a typo in equation (22).  That needs to be fixed before going ahead with the rest of the calculations.

 

@Neel What you have here is so far removed from what would be the right approach to the solution of the problem, that it leaves little to salvage.  In my earlier message I provided a general outline for solving the forced vibration problem.  What you have in your worksheet does not follow that outline at all.  The entire worksheet needs to be discarded and started over again.

But you should not feel bad about that.  My outline assumes some familiarity with orthogonal basis functions generated by solutions of boundary value problems, and how to expand an arbitrary function in such a basis.  It appears that you have not seen that material, so I wouldn't blame you for that.

I can provide the somewhat long and messy solution to the problem which you are attempting to solve, but without the proper mathematical background it won't be understandable.

Why are you interested in this problem?  If it is a student project, you should really begin with learning the basic mathematical ideas behind what is called eigenfunction expansions. Furthermore, the beam problem is one or two steps more advanced than a beginner problem in this subject, which would be:

First: The heat equation:
w_t= w_xx + f(x,t),
w(0,t)=0,   w(L,t)=0,
w(x,0) = g(x),

Second: The wave equation:
w_tt = w_xx + f(x,t),
w(0,t)=0,   w(L,t)=0,
w(x,0) = g(x),   w_t(x,0) = h(x).

You need to be able to confortably handle these before attempting the beam equation.  Do you have a textbook that shows how to solve these problems?  Learning that will be the first in your progression.  By the way, notice the initial conditions g(x), and h(x).  You will need to specify initial conditions to solve your beam problem as well.  I didn't see initial conditions in your worksheet.
 

@panke Specify the version of your Maple when you ask a question.  There is a pull-down menu just before you submit your question in which you enter the version.  Please don't ignore it.

The phase space of a 5th order ODE is five-dimensional.  The current version of Maple does three-dimensional pictures quite well.  You will have to wait for a future release of Maple to do plots in five dimensions.

@vercammen If you know about Maple's procs, then a better way would be through a proc, as in

make_prob := proc({n:=10, alpha:=-60,  beta:=20, P_val:=4})
    local Qs := alpha + beta*P;
    printf("There are %a identical firms.  "
           "Supply is given by Qs = %a.  "
           "When P = %a then quantity supplied is equal to %a.  "
           "Calculate the ...\n", n, Qs, P_val, eval(Qs, P=P_val));
end proc:

Then you may call the proc any number of times with the desired arguments.  Missing arguments take on default values as indicated in the first line of the proc's definition.  For example:

make_prob();   # all default args
There are 10 identical firms.  Supply is given by Qs = -60+20*P.  When P = 4 then quantity supplied is equal to 20.  Calculate the ...

make_prob(n=40, P_val=5);   # some optional args
There are 40 identical firms.  Supply is given by Qs = -60+20*P.  When P = 5 then quantity supplied is equal to 40.  Calculate the ...

 

@Neel There is no general way for that.  Each case needs to be treated on its own, but the steps are pretty much always the same..  So if you see what I did for the clamped-clamped case, you may do the same to analyze the clamped-free case (a cantilever beam).

Here is a pointer to simplify your calculations just a bit.

In what I presented earlier, I went along with general solution of the ODE produced by Maple as a linear combination of exponential and trigonometric functions with coefficients _C1 through _C4.  For our purposes that's not the best form of the general solution.  The first two terms may be replaced with hyperbolic functions, as in _C1*sinh(mu*x/L) + C_2*cosh(mu*x/L).  Since _C1 and _C2 are arbitrary, the changed expression is equivalent to the original.  I suggest that you make that change in the worksheet (manually) and then continue the calculations.  You will see that things work out more smoothly that way.

 

@Aminaple The idea suggested by Dr. Venkat Subramanian is quite good and I like it better than the solution that I presented earlier.  Here are the details of this alternative approach.

restart;

We convert the 2nd order differential equation to a system of 1st order differential
equations by introducing y(t) as the derivative of x(t)

de1 := diff(x(t),t)=y(t);

diff(x(t), t) = y(t)

Then replace "x ' (t) " and x*''*t with y(t) and "y ' (t)" in the differential equation and rearrange:

diff(y(t),t) + 2/25*y(t) + 4*x(t) = (1/25*diff(y(t),t) + 2/3*y(t) + 5/2*x(t)) * x(t-d):
isolate(%, diff(y(t),t)):
de_tmp1 := collect(%, x(t-d));

diff(y(t), t) = (((2/3)*y(t)+(5/2)*x(t))*x(t-d)-(2/25)*y(t)-4*x(t))/(1-(1/25)*x(t-d))

When t < d, the term x(t-d) picks up its values from the history, let call it "h(t)," so the differential

equation really has the form

de_tmp0 := subs(x(t-d)=h(t-d), de_tmp1);

diff(y(t), t) = (((2/3)*y(t)+(5/2)*x(t))*h(t-d)-(2/25)*y(t)-4*x(t))/(1-(1/25)*h(t-d))

Furthermore, when t < 0, we want y(t) to conform to the history, that is,  "y(t)=h ' (t)"  Therefore
"y ' (t) = h ' ' (t)."   We combine the three cases into a single differential equation:

de2 := diff(y(t),t) = piecewise(t<0, (D@@2)(h)(t), t<d, rhs(de_tmp0), rhs(de_tmp1));

de2 := diff(y(t), t) = piecewise(t < 0, ((D@@2)(h))(t), t < d, ((2*y(t)*(1/3)+5*x(t)*(1/2))*h(t-d)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*h(t-d)), ((2*y(t)*(1/3)+5*x(t)*(1/2))*x(t-d)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*x(t-d)))

The initial conditions at t = -d  are determined by the history function:

ic := x(-d) = h(-d), y(-d) = D(h)(-d);

x(-d) = h(-d), y(-d) = (D(h))(-d)

It's time to introduce the history function and the delay:

params := { h = cos,  d = 6/7 };

{d = 6/7, h = cos}

So our system looks like this:

sys := eval({de1, de2, ic}, params);

sys := {diff(x(t), t) = y(t), diff(y(t), t) = piecewise(t < 0, -cos(t), t < 6/7, ((2*y(t)*(1/3)+5*x(t)*(1/2))*cos(t-6/7)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*cos(t-6/7)), ((2*y(t)*(1/3)+5*x(t)*(1/2))*x(t-6/7)-2*y(t)*(1/25)-4*x(t))/(1-(1/25)*x(t-6/7))), x(-6/7) = cos(6/7), y(-6/7) = sin(6/7)}

We solve the system and plot the result

dsol := dsolve(sys, numeric):

plots:-odeplot(dsol, [t,x(t)], t=-eval(d, params)..12);

The graph agrees with what I had obtained with the earlier method.

 

Download mw2.mw

 

 

@Neel 

In the forced case you will have
"((&PartialD;)^2w)/((&PartialD;)^( )t^2)+1/(c^(2))((&PartialD;)^(4)w)/((&PartialD;)^( )x^(4))=f(x,t)".

To solve, express both w(x, t) and f(x, t) in terms of the modal functions

"w(x,t)=(&sum;)`alpha__j`(t)*`u__j`(x),  f(x,t)=(&sum;)`beta__j`(t)*`u__j`(x),"

where the coefficients `&beta;__j`(t) are known but `&alpha;__j`(t) are unknown.

Plug this into the PDE, and equate the coefficients of u__j(x).

That will give you an ODE of the form "`alpha__j`' ' + `h__j`*`alpha__j`=`beta__j`" for each j.
Solve those ODEs for `&alpha;__j`(t) and plug the result in the expression

w(x, t) given above.

 

 

 

@Neel What are the boundary condtions?  I cannot tell by reading through your worksheet.

@janhardo Your interpretation of my calculations is quite correct.  Well done!

@AHSAN Even after changing q(y)=j to u(y)=j the problem still makes no sense.  The "bc" and "bc1" in your formulation stand for "boundary condition" but you have not said where the boundaries are.

Go to where you found the problem, read it carefully, then come back with the correct statement.

 

I looked through your worksheets but was unable to understand them.  It will be very helpful to your readers to have some verbal description of the purpose of the worksheets.  You may want to supply the following information:

  1. It appears that you are solving the PDE 

    with the initial conditions u=(x,0)=sin(Pi*x), diff(u(x,t),t)=0 on the interval (0.1).  Is that correct?
  2. You state the boundary conditions as
    bc1:=u(x,t):        bc2:=u(x,t):
    I don't know that those mean.  Explain.
  3. It appears that you discretize the equation through centered finite differencing in the x direction. Is that correct?
  4. Then you discretize the equations in the time direction in some way.  That does not seem to be consistent with your stated objective of applying the method of lines, which approaches the problem through solving a system of ODEs, but perhaps I have misunderstood.  Explain.
  5. What is the idea behind your non-classical method?
  6. You write "I need to compare the result of the method with the exact solution".  What is the exact solution?
  7. You must know that you may solve your PDE with Maple's pdsolve() (numerically).  You are aware of that, aren't you?
  8. Why have you posted two worksheets?  What is the purpose of each?  Provide a guide.

 

@janhardo There are neither double integrals nor polar coordinates there.  The area is calculated by simple integration in Cartesian coordinates.

First 27 28 29 30 31 32 33 Last Page 29 of 91