## 6217 Reputation

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

• 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

## Not really...

The OP originally posted

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

simp38:=proc(f,x0,xn,nn)

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

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

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!

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

## Well theoretically...

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

de3.mw

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

## Things to consider...

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

## Fundamental point...

@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

## Makes little sense...

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

diff(u(x,t),t)+u(x,t)*diff(u(x,t),x)=0

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$

## Q cancels...

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:-(

## Agreed...

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

de2.mw

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