16 years, 0 days

## A symbolic solution...

It appears that Maple 2018 can find a symbolic solution. It is big, but integrating the first ode w.r.t. x and using ics helps and is done here:
Note added: The code below assumes that Q is a constant. If Q depends on x the method may still work. If not use numerical solution as Carl does.

```restart;
sys_ode:=    2*C__5*diff(w(x),x) - C__1*diff(u(x),x\$2) - 2*C__4*diff(w(x),x\$3) = Q,
2*C__2*u(x) - C__1*diff(w(x),x\$3) - 2*C__3*diff(u(x),x\$2) = 0;
ics:= w(0) = 0, D(u)(0) = 1, (D@@2)(w)(0) = 0, D(u)(100) = 1, (D@@2)(w)(100) = 2;
###
map(int,sys_ode[1],x=0..x,continuous);
ode1:=eval(%,{ics});
ode2:=sys_ode[2];
### Now solve symbolically:
res:=dsolve({ode1,ode2,ics}):
U,W:=op(subs(res,[u(x),w(x)])):
Digits:=200:
## Example plot:
plot(eval(U,{Q=1,C__1=1.2e8, C__2=2e9,C__3=2e8,C__4=0,C__5=3e7}),x=0..100,size=[2000,default]);

```

The plot of W is done similarly:

```plot(eval(W,{Q=1,C__1=1.2e8, C__2=2e9,C__3=2e8,C__4=0,C__5=3e7}),x=0..100,size=[2000,default]);
```

## Formal parameter...

In procB you have the formal parameter the_equation. That appears in the procedure body of procB in the print statement as well as in the call to procA. Wen procB is called with the equation y=3 the formal parameter the_equation will be replaced by y=3 in all its occurrences in the body of procB. Thus the call to procA will be procA('y=3' = (y=3)), so that the default NULL will be used by procA.

`procA(':-the_equation'=the_equation);`

## Inevitable...

Your sys1 will necessarily have a singularity for some t1 < tini.
The reason is that as long as a(t) > 0 it follows that a(t) is increasing. Going back in time a(t) will decrease to 0 at some finite value of t (called t1 above). Here a'(t1) = infinity.
To prove that that is what happens is rather straightforward. Your ode is of the form

`ode:=diff(a(t),t) = k*a(t)*sqrt(1/a(t)^3 + m); # k>0, m>0.`

with a(t0) = a0 > 0.

As long as a(t) > 0 we have

`diff(a(t),t) > k*a(t)*sqrt(1/a(t)^3) = k/sqrt(a(t))`

Rewrite this to get

`sqrt(a(t))*diff(a(t),t) > k`

Integrate from t to t0 (with t < t0) to get

`2/3*(a0^(3/2) - a(t)^(3/2)) > k*(t0-t)`

from which follows that

`t > t0 - a0^(3/2)/k*(2/3)`

since we are assuming that a(t) > 0. Thus a(t) will necessarily approach 0 from above at a finite value of t < t0.

## expand...

Try this where I have replaced the letter i by the letter I:

```restart;
A:=Matrix(6,(i,j)->cos(2*Pi/6*(i-1)*(j-1))-I*sin(2*Pi/6*(i-1)*(j-1)));
B:=conjugate~(A);
expand~(A.B);
```

## Eliminate one of the two...

If you eliminate one of CNO and CNH3 you can easily solve the system. From your two odes it follows that the difference between CNO and CNH3 is constant, thus in the case presented equal to zero.

```restart;
r := 2900*exp(-53300/(R*T))*CNO(t)^0.62*CNH3(t)^(-0.05);
ode := diff(CNO(t), t) = -r: #(negative because they are decaying with time)
ode2 := diff(CNH3(t), t) = -r:
ics := CNO(1020) = 1.6, CNH3(1020) = 1.6; #(sets up known initial conditions)
ode := subs(R = 8.314, T = 473, ode):
ode2 := subs(R = 8.314, T = 473, ode2):
#sys_ode := (ode, ode2) ;
odeS:=subs(CNH3(t)=CNO(t),ode);
#dsolve({odeS,CNO(1020)=1.6}); #??? Why does it take so long (or forever?)
sol:=dsolve(odeS);
eval[recurse](sol,{t=1020,CNO(1020)=1.6});
solve(%,{_C1});
eval(sol,%);
res:=solve(%,CNO(t)); # The solution for CNO(t)
t0:=solve(res=0,t);
plot(res,t=1020..t0);
```

The solution becomes 0 in finite time t0 (approximately 1775.4) due to the power of the right hand side of odeS being positive and less than 1. After that the solution continues as zero.

If the initial conditions are different it appears that you run in to an integral that Maple cannot find. A numerical solution may then be used:

```ics2 := CNO(1020) = 1.6, CNH3(1020) = 1; #Example of different initial conditions
# Now CNH3(t)-CN0(t) = constant = 1 - 1.6 = -0.6
odeS2:=subs(CNH3(t)=CNO(t)-0.6,ode);
res:=dsolve({odeS,CNO(1020)=1.6},numeric);
plots:-odeplot(res,[t,CNO(t)],1020..2000);
```

Notice though that since CNH3(t) = CNO(t) - 0.6 we must stop the integration when CNH3(t) = 0.
You could just plot CNH3(t) like this:

plots:-odeplot(res,[t,CNO(t)-0.6],1020..2000);

An alternative:

```resE:=dsolve({odeS,CNO(1020)=1.6},numeric,events=[[CNO(t)-0.6,halt]]);
resE(2000);```

You will be told that integration stopped at t = 1279.9381.

## Color...

Make the following addition and change:

```col:=["black", "blue", "red", "pink"];
for k to 4 do
R := dsolve(eval({bc, sys}, Ha = L[k]), fcns, type = numeric, AP);
AP := approxsoln = R;
p1u[k] := odeplot(R, [y, U(y)], 0 .. 1, numpoints = 100, labels = ["y", "U"], style = line, color = col[k])
end do:
```

## A way, but with problems...

Begin by revising the definition of F and then use a revised version of your second method.
## I revised my answer too  :-)
## Introduce Y'' as Yd2(x) (the easiest solution in this case, I think):

```F := dsolve({cond, sys,diff(Y(x),x,x)=Yd2(x)}, numeric, output = listprocedure);
Y1 := subs(F,Y(x));
Y2 := subs(F,diff(Y(x), x));
Y3 := subs(F,Yd2(x));

```

The remaining ones ( R(x), mu(x) , and mu'(x) ) are found quite similar to Y(x), and Y'(x).

## There appears to be singularity problems.

```F := dsolve({cond, sys}, numeric, output = listprocedure); #No Yd2
plots:-odeplot(F,[x,Y(x)],0..Pi);
```

You seem to be interested in x = 0..Pi, but you get from odeplot:

Warning, cannot evaluate the solution further left of 2.6226594, probably a singularity

Furthermore when trying to include diff(Y(x),x,x) = Yd2(x):

```F := dsolve({cond, sys,diff(Y(x),x,x)=Yd2(x)}, numeric, output = listprocedure);
plots:-odeplot(F,[x,Y(x)],0..Pi);
```

then we get an error I don't remember having ever seen before:

Error, (in dsolve/numeric/SC/IVPrun) invalid reversal of integration direction

Going the other way is OK at least until 3.66:

```plots:-odeplot(F,[x,Y(x)],Pi..3.66);
```

# When I use the second version of F (i.e. the one including Yd2) and try F(3) then I get this error:

Error, (in unknown) cannot evaluate the solution further left of 3.6605260, probably a singularity

which is surprising to me because who told F to go from Pi to 3.6605260 before going to 3?

## I notice that the weird error disappears if cot(x) is replaced by 1. Singularity problems don't disappear though.
###############
The problem (unrelated to other singularities) appears to have to do with the initial singularity at x = Pi where cot(x) is infinite (undefined):
Try

```F := dsolve({cond, sys}, numeric, output = listprocedure,abserr=1e-15,relerr=1e-13);
F(3);
```

and you get the error message:

Error, (in unknown) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

## Probably indeed a singularity!...

Actually, in Maple 2018 I get the complaint:

Warning, cannot evaluate the solution further right of 3.4774930, maxfun limit exceeded (see ?dsolve,maxfun for details)
So since maxfun can be set, I did this:

`res := dsolve(sys union ics, numeric, maxfun=10^6);`

But the complaint just changed the t-value slightly to 3.4775188.
There may be an actual singularity though.
Try this:

```sys:=solve(odes,diff~({x,y,z}(t),t)); # Solve for the derivatives
## Then do
plots:-odeplot(res, [t,rhs(sys[1])], Pi .. 4,caption="The derivative of x(t)");

```

What I see is this:

Using method=rosenbrock gives me the error you report:

`res := dsolve(sys union ics, numeric, maxfun=10^6,method=rosenbrock);`

## Function definition...

The correct way to define a function is like this:

```T:=x->1/2*x*sqrt(-x^2+(60-x)^2);
##
solve(diff(T(x),x)=0,x);
##
T(20);
```

## Calculus1:-ApproximateInt, NumericalAnal...

If you haven't already done so, you should take a look at Student:-Calculus1:-ApproximateInt and Student:-NumericalAnalysis:-Quadrature at least for comparison with your own code.

```restart;
fx:=x^2*sqrt(25-x^2);
Student:-Calculus1:-ApproximateInt(fx, x=0..5,method=trapezoid,partition=100);
evalf(%);
Student:-Calculus1:-ApproximateInt(fx, x=0..5,method=midpoint,partition=100);
evalf(%);
```

## Note: The following NumericalAnalysis:-Quadrature method is the same as the Calculus1:-ApproximateInt midpoint method:

```Student:-NumericalAnalysis:-Quadrature(fx,x=0..5,method=newtoncotes[open,0], partition=100);
```

## Corrections...

mmcdara already pointed out that you are missing the integration variable in the procedure pn.
But pn has another problem: pn has the formal parameter t  ( pn:=proc(i,n,t) ... ). OK.
But in the body of pn you have int(hd[i],t). That forces you to use literally t when using the procedure because hd[i] depends on the global t. This is bad (weird):

```hd:=t^2:
p:=proc(t)
int(hd,t)
end proc:
p(t);
p(s);
```

Simple solution for this case: leave out the third formal argument in pn.
##
You write:  " #alpha1(t) and alpha2(t) must be defined as functions of t ".
but you don't do that.
And what is f?
(I agree with mmcdara that you should use add instead of sum for concrete finite sums).
###########################
Putting in the functions from Example 2 in your reference just before defining RHS:

```f:=t->t^2;
alpha1:=t->-1:
alpha2:=t->t:
y0:=1:
y1:=-1:```

After that do for comparison with the numerical solution found by Maple (which is good):

```ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2;
res:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric);
plots:-odeplot(res,[t,Y(t)-y(t)],0..1);
```

The error in the method with the given choices:

With J=7 (so N=128) you get the error:

## Cannot reproduce...

When doing this:

```restart;
A:=LinearAlgebra:-RandomMatrix(5);
B:=LinearAlgebra:-RandomMatrix(5,datatype=float);
C:=LinearAlgebra:-RandomVector(50,datatype=float);
```

I see in the Variables list to the left:
A   5x5 Matrix
B   5X5 Matrix
C   50-element Vector[column]

###

`E:=Array(1..100,i->i);`

E   100-element Array.
###
`Standard Worksheet Interface, Maple 2018.0, Windows 10, March 10 2018 Build ID 1298750`

## convert/Vector...

Your F is a Matrix with one column, so you need to convert it to a Vector since the left hand side of sys1 is a Vector.

```whattype(F); # Matrix
op(1,F); # 4 rows, 1 column
##
sys1:= diff~(XX,t)+A.XX=convert(F,Vector):
dsolve(sys1); # returns a DESol structure
```

The DESol structure is basically just a rewriting of the input.
Your best bet is numerical solution, but then you cannot input a vector equation. Rewrite as a system using Equate:

`sys:=convert(Equate(lhs(sys1),rhs(sys1)),set);`

For numerical solution you need concrete values for all parameters  a_0, i, m_0, t, v, v_0 (and initial values).

## System of odes...

I assume that the equations you are talking about are f=0, g=0, h=0.
If so you can turn those into odes:

```restart;
eqs:={t^2*x^2 + t*x + 2*x - 1=0, t^2*y^3 + t*y^2 + 2*y - 1=0, 2*t^2*z^3 + t*z + 3*z - t^2=0};
eqs2:=subs(x=x(t),y=y(t),z=z(t),eqs);
eval(eqs2,t=0);
ics:=solve(%,{x,y,z}(0));
odes:=diff~(eqs2,t);
res:=dsolve(odes union ics, numeric);
plots:-odeplot(res,[x(t),y(t),z(t)],0..4); p1:=%:
```

There will be at least one other curve in this case if you allow initial points for other t values as well. Using t=1 gives two solutions (one in addition to the one we found above):

```eval(eqs,t=1);
ics1:=solve(%,{x,y,z});
evalf(allvalues(ics1));
remove(has,[%],I);
ics2:=subs(x=x(1),y=y(1),z=z(1),%[2]);
res2:=dsolve(odes union ics2, numeric);
plots:-odeplot(res2,[x(t),y(t),z(t)],0.2..4); p2:=%:
plots:-display(p1,p2);
```

There may be a smarter method.

## evalc...

evalc(Re(expr)) will find the real part of expr. evalc will assume that all variables represent real numbers.
Thus you can do:

```restart;
lambda2 := a/2-3/2+1/2*sqrt(a^2-10*a-23);
solve(lambda2<=0,{a});
solve(evalc(Re(lambda2))<=0,{a});
plots:-complexplot(lambda2,a=-8..3,thickness=3);
```

You can do the same with lambda3 := a/2-3/2-1/2*sqrt(a^2-10*a-23) of course.

 First 10 11 12 13 14 15 16 Last Page 12 of 149
﻿