Dr. Patrick T

## 2083 Reputation

13 years, 259 days

## maybe assuming and AllSolutions...

making assumptions about x may eliminate some candidate solutions, you may also  remove rootofs with the AllSolutions option in solve, maybe something like

```E[i, j] := Int(f[i](t)*u[j](t), t = -infinity .. infinity);
eqn[i, j] := u[j](x) = E[i, j];
k[i, j] := solve(eqn[i, j], x, AllSolutions) assuming x::positive;
```

## subs...

also does the trick (but couldn't tell you which if any is to be commended)

```subs(Result,alpha);
-1.536774530
```

## manual method...

There must be a way to automate the following (with the use of the option discont=true), and someone else will certainly show how:

```restart;
with(plots):
with(plottools):

fdiscont((x^2-4*x+5)/(x-2),x=1..3);

[1.99972287449048180 .. 2.00074755004038485]

```

fdiscont helps you find the discontinuities of your function -- okay, so the asymptote is near 2 -- weird that Maple doesn't understand it's exactly 2.

Then I split the plot into two parts, so that Maple won't show a vertical line at the asymptote -- I then create my own asymptote by hand and give it the look I want (including dashes):

```p1:=plot((x^2-4*x+5)/(x-2), x=1..1.99):
p2:=plot((x^2-4*x+5)/(x-2), x=2.001..3):
p3:=curve([[2,-10], [2,10]], color=red, linestyle=dash, thickness=2):
display({p1,p2,p3}, view=[1..3,-10..10]);
```

## maybe...

assuming you meant to type * instead of = in the formula for H:

```H:= (cos(n*L)*cosh(n*L))/(sin(n*L)-sinh(n*L)):
phi:=x->sin(n*x)-sinh(n*x)+H*(cos(n*x)-cosh(n*x)):
int(phi(x), x=0..L);
1/2 (4 exp(n L) sin(n L) - 4 exp(n L) sinh(n L)

- 2 cos(n L) sin(n L) exp(n L)

- cos(n L) cosh(n L) exp(2 n L)

+ 2 cos(n L) sinh(n L) exp(n L) + cos(n L) cosh(n L)

- exp(2 n L) sin(n L) + exp(2 n L) sinh(n L) - sin(n L)

+ sinh(n L) + 2 cos(n L) cosh(n L) sin(n L) exp(n L))

exp(-n L)/((sin(n L) - sinh(n L)) n)

```

## exponential notation on the vertical axi...

loglogplot(10^x, x=1..2);

on my Windows/Maple13/Classic system I do get exponentials on the vertical axis.

You did not give us your code.

## not very likely...

Explicit solutions are the exception rather than the rule. The following can be used to extract information from Maple's dsolve command.

```restart;
with(DEtools);
infolevel[dsolve] := 4;
```
```dsys1 := {diff(u(t), t) = -(3/2-sqrt(2))*u(t)+5*v(t), diff(v(t), t) = -5*u(t)-(3/2+sqrt(2))*v(t)-20*cos(10*t)*w(t), diff(w(t), t) = -.5-3*w(t)+20*cos(10*t)*v(t)};

dsolve(dsys1);

-> Solving each unknown as a function of the next ones using the order: [u(t), w(t), v(t)]
-> Calling odsolve with the ODE, diff(diff(diff(y(x),x),x),x) = -1/4*(-48000*cos(10*x)^2*y(x)*sin(10*x)+7980*cos(10*x)^2+3100*diff(y(x),x)*sin(10*x)-728*diff(diff(y(x),x),x)*cos(10*x)+600*diff(diff(y(x),x),x)*sin(10*x)-319200*cos(10*x)^3*y(x)-19291*y(x)*cos(10*x)+4800*cos(10*x)^3*diff(y(x),x)+9090*y(x)*sin(10*x)+32000*cos(10*x)^2*diff(y(x),x)*sin(10*x)-960000*cos(10*x)*y(x)*sin(10*x)^2-1989*diff(y(x),x)*cos(10*x)+16000*sin(10*x)^2+(1600*cos(10*x)*sin(10*x)-96000*cos(10*x)^2*y(x)*sin(10*x)+80*diff(diff(y(x),x),x)*sin(10*x)+240*diff(y(x),x)*sin(10*x)+2020*y(x)*sin(10*x)+3200*cos(10*x)^3*diff(y(x),x)+48*diff(diff(y(x),x),x)*cos(10*x)+274*diff(y(x),x)*cos(10*x)+606*y(x)*cos(10*x))*2^(1/2))/(3*cos(10*x)+20*sin(10*x)+2*2^(1/2)*cos(10*x)), y(x), singsol = none
Methods for third order ODEs:
--- Trying classification methods ---
trying high order exact linear fully integrable
trying differential order: 3; linear nonhomogeneous with symmetry [0,1]
trying high order linear exact nonhomogeneous
trying differential order: 3; missing the dependent variable
-> pFq: Equivalence to the 3F2 or one of its 3 confluent cases under 'a power @ Moebius'
trying a solution in terms of MeijerG functions
-> Try computing a Rational Normal Form for the given ODE...
<- unable to resolve the Equivalence to a Rational Normal Form
checking substitution methods
trying differential order: 3; missing the dependent variable
-> pFq: Equivalence to the 3F2 or one of its 3 confluent cases under 'a power @ Moebius'
trying a solution in terms of MeijerG functions
-> Try computing a Rational Normal Form for the given ODE...
<- unable to resolve the Equivalence to a Rational Normal Form
checking substitution methods
trying differential order: 3; missing the dependent variable
checking if the LODE is of Euler type
trying Louvillian solutions for 3rd order ODEs, imprimitive case
-> pFq: Equivalence to the 3F2 or one of its 3 confluent cases under 'a power @ Moebius'```

Maple then gives you an equivalent representation of your system in terms of one third-degree ODE in v(t). Run the above code and you can see that it's unwieldy.

I don't get the memory allocation error problem.

I am guessing that EU must be something like: Px*U(x)+Py*U(y), where U() is your utility function.

if your utility function is quadratic, since your constraint linear, you should/must be able to solve the problem in closed form.

if U(x)=1-x^2, say...

## Axel's contourplot with numpoints=50^3...

Takes a short while to produce but looks pretty smooth:

```contourplot( max(min(-z*(1/(x^2+z^2)^(3/2)-1), 1000),-1000),
x=-0.3..0.3, z=-0.3..0.3, numpoints=50^3,axes=boxed);
``` ## explicit use of max and min...

The max min trick is pretty cool. Its usefulness is particularly apparent, in this example, for small values of numpoints:

```plot3d( max(min(-z*(1/(x^2+z^2)^(3/2)-1), 1000),-1000),x=-0.3..0.3,z=-0.3..0.3,
numpoints=50,axes=boxed,lightmodel=light2,orientation=[27, 54]);
```
``` ```

versus

```plot3d( -z*(1/(x^2+z^2)^(3/2)-1), x=-0.3..0.3, z=-0.3..0.3,
view=[DEFAULT,DEFAULT,-1000..1000], numpoints=50,
axes=boxed,lightmodel=light2,orientation=[27, 54]); ```
```
```

Okay, I confess, I'm really practicing my image-posting skill.

## oops, sorry!...

you're absolutely right, I copied the wrong line of code, apologies.

## system explodes...

I think. You can see that by gradually raising the range of t from, say, 0..1 to 0..4, the following sums it up:

```plots[odeplot](sol,[t,Theta1(t)],0..2);
plots[odeplot](sol,[t,Theta2(t)],0..2);
```

## solve...

in this case, solve does it:

y=(-4+4*d)*x^2+(4-4*d)*x-3*d+2;
solve(%, x);

## numpoints=n, n>50...

one way to do it: numpoints

the default numpoints is 50, increase it for better effect:

contourplot(-z*(1/(x^2+z^2)^(3/2)-1),x=-3..3,z=-3..3,grid=[50,50],contours=10,numpoints=100);

## is this what you're looking for?...

with(plots):
f := x -> x^(-2):
logplot(f(x), x=0..2);

﻿