6217 Reputation

17 Badges

10 years, 99 days

MaplePrimes Activity

These are replies submitted by tomleslie

Agree with Dr. Venkat Subramanian, but adding a couple of further issues you have to address

  • The statement i <> j;  does nothing! It will return true or false, but since nothing subsequently depends on the truth or falsity of this statement, then it is redundant. Maybe(?) you mean something like

if i <> j;
then d1 := diff(x(z), z) = -G*x(z)*y[i](z)/IC-alpha*x(z);
       d2 := diff(y[i](z), z) = G*y[i](z)*y[j](z)/IC-alpha*y[i](z);
       dsys := {d1, d2};
       F := dsolve({ICon, op(dsys)}, [x(z), y[i](z)], numeric);

  • Since you are computing 2500 solutions (well maybe 2450 - see the first bullet) for 2500 (2450) different(?) differential equations you should be aware that as written each value of dsolve({ICon, op(dsys)}, [x(z), y[i](z)], numeric); will overwrite the previous one and at the end of the exercise, F will contain the solution to one differential equation with i=50 and j=50 (or possibly i=50 and j=49 - see  first bullet). Seems like an awful lot of work just to throw away the answers


The OP originally posted

 \int_{a}^{b} f(x) \, dx \approx \tfrac{3h}{8}\left[f(x_0) + 3f(x_1) + 3f(x_2) + 2f(x_3) + 3f(x_4) + 3f(x_5) + 2f(x_6) + ... + f(x_n)\right] .

I want to write a procedure for the simpson 3/8 rule using the above formula I took from wikipedia :

Unfortunately (s)he missed the next line in the Wikipedia entry which states that (emphasis added)

Note, we can only use this if n is a multiple of three.

The OP's code does not enforce this requirement (so neither does yours). Since either code accepts any number of points, neither is implementing Simpson's 3/8 rule, unless the number of points happens to be a multiple of 3

A trivial amendment to either code will ensure that it will always implement Simpson's 3/8 rule: a modified version of your code is shown below


local  n:=3*nn, h:=(xn-x0)/n,  x:=i->x0+i*h;

3*h/8*(f(x0)+f(xn)+3*add(f(x(1+3*i)), i=0..(n-2)/3)+3*add(f(x(2+3*i)), i=0..(n-3)/3)+2*add(f(x(3+3*i)), i=0..(n-4)/3));

end proc:

A similar modification will also work for the OP's code. With this modification, one can compare with maple's built-in command for an arbitrary function with arbitrary end-points using

kit:=n-> simplify(simp38(f,a,b,n));
maple:=n->Student[Calculus1][ApproximateInt](f(x), a..b, method=simpson[3/8], partition=n);
seq( testeq(kit(j)=maple(j)), j=1..10)

Part(a) can be proved by induction on k although I think the limits on the sum should be 1,k. Consider what the OP's original gives when k=1!

   A:=k->sum( Q*exp(-n*c*T), n=1..k):
# Check base case for induction
# test inductive step
   testeq( A(k)*exp(-c*T)+A(1) = A(k+1) );


If you take your 2 differential equations (or do all 5 if you want), and define odesys:={de1,de2}, then use odetest(puta, odesys), you will get a list of 2 (or 5) residuals. You can then use Preben's technique (solve/identity) to generate sets of solutions for each of the residuals.

What you are then looking for is an element which is common to both (or all 5) sets. For a variety of reasons this might be a bit ugly to do (possible, but ugly). However while experimenting with this, I discovered that one of your residuals is only ever satisfied when lambda=0, and the other would seem to have only non-zero values (unless one can require a value for w - ie treat is as a parameter in the same way as a,alpha, beta,gamma1,lambda).

See attached worksheet

So either

  1. typo(s) in ODEs
  2. no solution with given trial functions
  3. 'w' has to be treated as a parameter in residuals

You have three second-order differential equations. As a general rule of thumb this would require 6 initial conditions, and you only have 4 - not good :-(

Your initial conditions contain functions such as u(R(z)/R[0]) = 0: what is R(z)/R[0]??? Possibly R(z)/R(0), note parentheses: but then what does R(z)/R(0) mean?? basically that u(some_unknown_value)=0. As an "initial condition", that is about as useful as a chocolate teapot

Your differential equations contain the variable 'i' in a few places. Whilst 'i' is perfectly acceptable as a variable name, I hope you don't mean the square root of -1, which in MapleSpeak would be I


@Preben Alsholm

As Preben has pointed out a couple of posts ago, the objective of the odetest()  command is to verify whether a given list/set of solutions are valid for the given list/set of differential equations - check the manual which says

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.

The fact that your first application of the odetest command, ie

odetest(puta, odesys)

does not return 0, indicates that your putative solutions (puta) are not valid for your ODE, odesys.

I think you have to resolve this issue before you can do anything else. Either there is a problem with your ODE (typo, maybe), or your guessed(?) solutions are not in fact correct


Sorry but your question doesn't make a lot of sense to me. I don't really  understand the first part (problem set-up) at all, and fail to see how this relates to the second part (all the implicitplots). A few simple questions

When you say $u_[t] + uu_[x] = 0$, I assume the '$' symbols are completely meaningless. The '$' construct can do useful things in Maple, but I am reasonably sure that this is not what you want, so you actually mean u_[t] + uu_[x] = 0. Now I am pretty sure that you don't want indexed values (ie [] brackets), but I am also pretty sure that you don't want two simple functions u_(), and uu_() so I am going to make another wild guess that this is actually your differential equation, so what you mean by this "expression" would actually be represented in Maple as


which, if memory serves, looks a bit like the Hopf equation often written as ut+u*ux=0 (althoughMaple does not allow the use of subscripts to indicate partial differentiation - you can "spoof" it but I wouldn't bother)

Now I have to work out what you mean by

$u(x,t) = a$ if $x < -1$
           = $b$ if $ -1 < x < 1$
            = $c$ if $ x > 1$

Again I'm pretty sure that you don't need/want all of those '$' symbols - so the first step is to produce

u(x,t) = a if x < -1
         = b if  -1 < x < 1
         = c if  x > 1

Defining a function like this in Maple is best done by using the piecewise construct

piecewise( x<-1,                 u(x,t)=a,
                x>=-1 and x<1, u(x,t)=b,
                x>=1,                u(x,t)=c
I'm now guessing that in fact this expression is intended to determine your initial conditions (since you claim to have an initial value problem) - but the above does not define any "initial" values, so my next guess is that you mean to define the initial conditions

ics:= piecewise( x<-1,                 u(x,0)=a,
                       x>=-1 and x<1,  u(x,0)=b,
                       x>=1,                u(x,0)=c

So the summary of the problem would be that you want to solve the partial differential equation

diff(u(x,t),t)+u(x,t)*diff(u(x,t),x)=0 (whhihc I think is the Hopf equation)

subject to the initial conditions

ics:= piecewise( x<-1,                 u(x,0)=a,
                       x>=-1 and x<1,  u(x,0)=b,
                       x>=1,                u(x,0)=c

Before I go any further (if I go any further), please confirm whether the above is correct, or if I have been wasting my time piling wild guess on wild guess.


I have no idea what you are trying to do with those implicitplot commands or how they relate to the first part of the problem - I would assume that you are trying to plot the solutions to the pde but can make absolutely no sense of what is written here


last time, when I said

"I suggest you post a worksheet showing your calculation so far with a detailed description of which part is failing"

I meant it. You have to help us, to help you

Don't copy-and-paste. Use the big green up-arrow (right-hand element in the second row of icons in the toolbar) to upload your worksheet

As a general rule I don't do homework problems unless the original poster can show some signs of having at least tried to solve the problem themselves. Even if I felt inclined to solve this problem, I could not so (neither could anyone else) given the supplied information, so you will have to try harder.

The MOL technique relies on discretizing one all but one dimension: the reaming dimension is treated as continuous. From your comment

t=1,delta t=0.1,x=1,Nx=10

which of x and t are considered discrete and which is continuous?? Does t go from 0 to 1 in steps of 0.1?? This is OK but unusual because in most MOL problems t is considered continuous and spatial dimensions are discretized. You can do it the other way round, , but I would want some definite confirmation before I even started. Other wise I could spend time solving the wrong problem - so which is it??

In your definition of the differential equations you have

du/dt=0.01 d^2u/dx^2

Does d^2u/dx^2 equal anything in particular - you know like 0, 1, exp(x)*cos(t)+myLatestShoppingBill?? Problem is insoluble without this information

Like I said - Try Harder

@Alejandro Jakubi 

After no response on this site, I submitted it to Maple Technical support and the last I heard from them (26Aug2014) was

Hello Tom,
Thank you for contacting Maplesoft. I suspect this behaviour may be due to the thread safe ability of the seq command.
I have forwarded this problem to our second level support for review. I will keep you informed as to the outcome.
Technical Support Analyst

so I never did find out why $ was so much slower than seq - I just stopped using $

@Preben Alsholm 

Preben is completely correct. I got so bogged down fixing the simple details of you problem that I lost sight of your ultimate objective - mea culpa

When using the odetest() command, you are supplying a list(set) of differential equations and a list(set) of functions, and asking Maple to verify whether the functions actually are *solutions* of the differential equations. If they are, then odetest will return {0=0} (or maybe[0=0]).

So with the amendments and comments which have been made so far, progress may have been made, but the supplied function list is a long way from being a solution of your differential equation.

The only combination of the remaining parameters I have come uo with so far whihc does form a solution is {a=0, gamma=beta}. Since this is equivalent to making A, B, C. phi independent of r, I have a funny feeling this isn't a solution you want!!


@Rouben Rostamian  

I agree with you: however I actually remember the very first post I ever made on this site which was to ask why repeated calculations using the $ construct ran about 10X slower than those using the seq construct. I never did get a satisfactory answer Maple techies eventually decided it had something to do with ensuring thread safety: I didn't particularly buy this (still don't), but since I discovered the time penaly involved, I have prety much given up on the $ construct


Sorry but Q does seem to cancel.

You can see why if you remove phi(r) from the substitution list: ie use

puta := {A(r) = (1-a/r)^alpha, B(r) = (1-a/r)^beta, C(r) = (1-a/r)^gamma};

followed by

collect(collect(odetest(puta, odesys), phi(r)), diff(phi(r), r));

which will "gather" terms in phi(r) and diff(phi(r),r).

There are only two terms containing these functions. The first depends on (diff( phi(r), r) )^2/phi(r)^2, and the second depends on diff(phi(r),r)/phi(r). The premultiplying factor of Q in the expression for phi(r) will give Q^2/Q^2 in the first term and Q/Q in the second one - so Q always cancels.

If you *know* that your system should depend on Q then I can only conclude that you have typo in your original equation: can't help with that:-(

@Preben Alsholm 

Attached is a corrected worksheet with all functions correctly defined. As well as, for example, A when you mean A(r), there are also a few places where you have A (r) (note the space) which also need correction.

The last couple of commands in the attatched worksheet show a couple of quick attempts to "collect" powers of r or (1-a/r)^.., nether of them particularly successful at "prettifying" the original expression

First 145 146 147 148 149 150 151 Last Page 147 of 156