Preben Alsholm

13728 Reputation

22 Badges

20 years, 241 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

It seems to turn up when using 2d-input. When using 1d-input ("Maple Notation") it doesn't show up.

Since the other replies use numerical solution of the ode, I thought I would point out that of course you can plot the exact solutions, when they are available, which admittedly often is not the case. But certainly in the example you mention. To make the plots more exciting I have used animation.

with(plots):
de := diff(y(x),x,x) + 4*y(x) = sin(2*x);
res:=dsolve({de,y(0)=y0,D(y)(0)=y1});
animate(plot,[eval(rhs(res),y0=0),x=0..3*Pi,caption=typeset([y(0)=0,D(y)(0)=y1])],y1=-2..2,trace=24,paraminfo=false);

#And now an animation inside another animation:
animate(animate,[plot,[rhs(res),x=0..3*Pi,caption=typeset([y(0)=y0,D(y)(0)=y1])],y0=-2..2,frames=10,paraminfo=false],y1=-2..2,frames=10,paraminfo=false);

I'm pretty sure that the problem is your use of !!! (Execute the whole worksheet).

Remember that it is the last assignment (timewise) that counts. If after executing the whole worksheet you return to a line above, then Maple uses that last assignment which may not correspond to the situation in that line.

Executing the whole worksheet is very often not a good idea, unless you have no intention to return to a previous line. 

If you do

dsolve(ode);

Maple will find the same solutions you could find by using the usual trick of multiplying the ode by y' and integrating to obtain

(*) diff(y(x),x) = (+/-) sqrt(a*y(x)^2-1/2*y(x)^4+1/3*y(x)^6+C)

Using the boundary condition D(y)(0) = 0 gives you a connection between C and y(0). Using y(2) = 0 when separating variables in (*) gives you an equation between y(0) and a by integrating from x = 0 to x = 2, something like

int(1/sqrt(a*y^2-y^4/2+y^6/3-(a*y0^2-y0^4/2+y0^6/3)),y=y0..0)=2; (using the plus sign from (*)).

Now you have to do integrate (!!!) and solve for y0 in terms of a. That is most likely not possible.

Even solving for y0 numerically when given a seems unlikely to succeed.

Do you have any reason to believe that there are nontrivial solutions to your boundary value problem for some value of a (or all)?

Looking into the numerical solutions it seems that there are solutions for negative values of a, e.g. taking a = -0.5 then there are solutions for y0 near -0.43 and -1.

Doing

dsolve({eval(ode,a=-0.5),D(y)(0) = 0, y(0) = -1},numeric);
plots:-odeplot(%,[x,y(x)],0..2);

seems to confirm that.

To answer your last question:

Solving D(L)(z) =( L(z)/(1+z) )+( (1+z)/H ), with H=Y(z) being a known numerically determined function, you can use the option known = Y.

You first have to make sure that Y returns unevaluated when not receiving numeric input (that is  'procname(z)' below).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=1) end if:
p := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)/Y(z),L(0)=0}, type=numeric, range=0..5,known=Y);
odeplot(p,[z,L(z)],0..5);

It might be better or at least as good to take the previously mentioned differential equations approach and now solving a system of two ode's numerically:

yp := implicitdiff(eq(z), H, z);
ode:=diff(H(z),z)=subs(H=H(z),yp);
p := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)/H(z), ode,L(0)=0,H(0)=1}, type=numeric, range=0..5);
odeplot(p,[z,L(z)],0..5);

I looked at your worksheet, but didn't run it.

However, just before your problems start you define C to be IdentityMatrix(2). This is very likely to cause all kinds of problems because C has already been used as C[alpha.f].

It is never a good idea to use a variable name and an indexed version of it independently.

Try this:

restart;
allvars:={C[alpha.f]=20,C[alpha.r]=24}:
C:=LinearAlgebra:-IdentityMatrix(2):
allvars;
Error, bad index into Matrix

Why start another thread and make it difficult for people to follow?

Anyway, the syntax is

implicitplot(eq(x), x=0..10,y=0..120);

and gives a curve above the x-axis. I am assuming that, as previously,

eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;

Using the example from one of your other questions you can use either fdiff or implicitdiff to plot anything with diff(y(x),x) in it.

To plot y(x) + 1 is straightforward.

eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;
Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);
plot(Y+1,-10..10);
plot(fdiff(Y,[1],x),x=-2..2);
yp:=implicitdiff(eq(x),y,x);
plot(eval(yp,y='Y(x)'),x=-2..2);

Rewriting for convenience:

eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;

The numerical solution procedure using fsolve:

Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);

Did you try fsolve?

fsolve({fff},{T, x[22], x[23], x[31], x[32], y[21]});

With a starting point (maybe you expect T > 0?)

fsolve({fff},{T=12, x[22], x[23], x[31], x[32], y[21]});

u1:=sin(sqrt(K)*t/sqrt(N));
subs(K=N*omega^2,u1);
simplify(%) assuming positive;

#No problem with the next version

u2:=sin(sqrt(K/N)*t);
subs(sqrt(K/N)=omega,u2);

 

The necessity of assumptions is due to the fact that sqrt(K/N) is not always equal to sqrt(K)/sqrt(N).

As an example consider K=1 and N = -1, in which case sqrt(K/N) evaluates to i, whereas sqrt(K)/sqrt(N) evaluates to 1/i = -i.

I guess instead of d = 1 you meant degree=1?

curvefit := Spline(roi, x, degree = 1):

#Then solve works:

solve(curvefit = .3, x);

# Since curvefit is not a polynomial only one solution is returned by fsolve at a time:

fsolve(curvefit = .3, x);


# but you can add a starting point to find the other solution:

fsolve(curvefit = .3, x = 740);

#You could also try
Student:-Calculus1:-Roots(curvefit = .3, x = 685 .. 761);

LinearAlgebra:-Norm(<-1,3>-<4,-2>,2);

I am assuming that you don't want to just plot the solutions to the equation, but rather e.g. the left and right hand sides:

eq:=x^4+x^2=3*x^2+8;
solve(eq,x);
plot([lhs(eq),rhs(eq)],x=-3..3,0..30);

#or you could plot the difference:

plot(lhs(eq)-rhs(eq),x=-3..3);
#Alternatively:
plot((lhs-rhs)(eq),x=-3..3);

You could write the points as lists (using square brackets [ ]) or as column vectors (using < >):

([0,2]+[1,3])/2;
(<0,2>+<1,3>)/2;
([a,b]+[c,d])/2;
(<a,b>+<c,d>)/2;

First 147 148 149 150 151 152 153 Last Page 149 of 160