11. Solving Differential Equations  (and Systems)

2015-05-29

$Version

11.1  Analytical Solution

K. Rektorys:  Survey of applicable Mathematics.  MIT Press, 1969 contains a long list of ordinary differential equations with solutions. 
E. Kamke, Differentialgleichungen. Lösungsmethoden u. Lösungen. [Bd] 1.2 Stuttgart: Teubner 1979-1983, 
	list ordinary and partial differential equations together with their solutions. 
A.D. Polyanin, V.F. Zaitsev, Handbuch der linearen Differentialgleichungen.	
Spektrum Akademischer Verlag, 1996.

For many but not all ordinary differential equations enumerated in these books analytical solutions can be obtained by the following commands.

11.1.1 Solving a single differential equation

DSolve[eqn, y[x], x] 	render the general solution of the differential equation given in eqn, i.e.   			 ---------		y[x],  taking  x  as the independent variable
					 						              	       		 
DSolve[{eqn,indat},y[x], x] render the solution belonging to initial  data given in  indat.

General solution:  (C[1], C[2], C[3], ...  are  the integration constants.)

11.1.1.1	  Some examples

Clear[x, y, a]
DSolve[ y'[x] == a y[x], y[x], x]

DSolve[ y''[x] + k^2 y[x] == 0, y[x], x]

DSolve[ y''''[x] + k^4 y[x] == 0, y[x], x]

DSolve[ y''[x] + y'[x]/x  +  (1 - n^2/x^2) y[x] == 0, 
		y[x], x ]

This is Bessel's differential equation. If  n  is not an integer, so n = ν, the two functions  (Subscript[J, ν])(x) and (Subscript[J, -ν])(x)
are linealy independent. Howerver, for  n = integer, these two solutions are dependent; so the Neumann 
solution  (Subscript[Y, n])(x) must be used as a second solution.

DSolve[ y''[x] + y'[x]/x  +  (1 - 4/x^2) y[x] == 0, 
		y[x], x ]

DSolve[ y''[x] + 2 y'[x]/x + (1 - n(n + 1)/x^2) y[x]==0, 
		y[x],x ]//ExpandAll

The above differential  equation  is obtained  when the Hemholtz equation  (Δ + (k^2))ψ = 0 is solved  in spherical  coordinates  by separating  the variables   ψ(r,θ,ϕ) := R(r ) Θ(θ) Φ(ϕ).  The resulting  equation for the radial  variable  R(r)  is the above differential equation with  (k^2) in place of 1. So the radial solution  is R(r)  = (Subscript[z, n])(k r), where  the spherical Bessel  function,  (Subscript[z, n])(x) =  (Sqrt[2/πx ])(Subscript[Z, n +1/2])(x), 
 is propotional to a common Bessel function of order  n + 1/2.

so = DSolve[ y''''[x] +3 y'''[x] +  k^4 y[x] == 0, y[x], x]

A linear differential  equation with constant coefficients  as given above can be solved by an  exponential  (e^( α x)),
where  α must be a  root of  the characteristic   polynomia, which in the  present case is   (α^4)  + 3 (α^3) + (k^4) .

In quantum mechanics the eigenvalue equation of the one-dimensional harmonic oscillator is given by the following second order differential equation, see also § 11.2.7.

seh = DSolve[ y''[x] + (ε - x^2) y[x] == 0, y[x], x]

sh = DSolve[ y''[x] - 2 x y'[x]+ 2 n y[x] == 0, y[x], x]

For  nonnegative integers  n, the function  HermiteH[n,x] beomes a  polynomial and the  product 
(e^(-(x^2)/2))  HermiteH[n,x]  fulfils the  boundary conditions   y(±∞) = 0.   This fixes the  dimensionless
 energy parameter  ε  by  ( ε -1) / 2  = 2n , n = 1,3,5, ...  ; i.e.  ε = k + 1/2, k = 0, 1, 2, ... .

11.1.1.2	  Solution Matching Given Data as eg. in Initial Value Problems

DSolve[ {y'[x] == a y[x], y[0] == y0}, y[x], x]

DSolve[{y'[x] - a y[x] == d Exp[c x], y[0] == y0}, y[x], x]

DSolve[{x[0]==3, x'[0]==0, x''[t]+x'[t]-x[t]==6},x[t],t]//Flatten

The data need not be initial data;  for example,  the data may be the value of the wanted solution at two places:

DSolve[{x[0]==3, x[5]==9, x''[t]+x'[t]-x[t]==6},x[t],t]//Flatten

11.1.1.3	  Plotting Solutions

Clear[x,y,f];
so = DSolve[ {y'[x] == 2 y[x], y[0] == 3}, y[x], x]//Flatten

g = y[x] /. so

Plot[g,{x,0,4},  AxesLabel -> {"x", "y[x] = 3 \!\(\*SuperscriptBox[\(E\), \(2\\ x\)]\)"}, ImageSize -> 200]

11.1.1.4   Printing  Numerical Values of Solutions

tso = Table[{x,g}//Flatten,{x,0,4,.5}]

TableForm[{"x   y(x)"}]
TableForm[tso]

11.1.2  Solving for solutions expressed as pure functions

DSolve[eqn,  y  , x]  	 give a solution for y in pure function form
		      ---

Note the difference in the calls of NDSolve[] above and at the beginning of  11.1.1.  This has farreaching consequences for the storage of the result and the reaction on further commands. When the solution obtained as a common function is inserted in an expression requiring the solution, the substitution is only done for the function, but not for the derivative nor for the inital values. Whereas for a  solution obtained as a pure function, the substitution is done for the function, the derivative(s) and the intial data. The difference is shown below.

Clear[x,y,f];
f = DSolve[ {y'[x] == 2 y[x], y[0] == 3}, y[x], x]//Flatten

y[x] + 2 y'[x] + y[0] /. f

Clear[x,y,f];
f = DSolve[ { y'[x] == a y[x] }, y, x]//Flatten

y[x] + 2 y'[x] + y[0] /. f

Clear[x,y,h];
h = DSolve[ {y'[x] == 2 y[x], y[0] == 3}, y, x]//Flatten

y[x] + 2 y'[x] + y[0] /. h

g = f /. {a -> 2, C[1] -> 3}

Plot[y[x]/.g,{x,0,4},  AxesLabel -> {"x", "y[x] = 3 \!\(\*SuperscriptBox[\(E\), \(2\\ x\)]\)"}, ImageSize -> 200];

gives same figure as above.

Clear[y,x]
syseqn = { y''[x] + k^2 y[x] == 0}

DSolve[syseqn, y[x], x]

Clear[y,x,f,g]
syseqn = { y''[x] + k^2 y[x] == 0, y[0] == 0, y'[0] == 1 }

so = DSolve[syseqn, y[x], x] //Flatten

su = DSolve[syseqn, y, x]//Flatten

ExpandAll[su]

g = y[x] /. su

11.1.3   Solving Systems of Differential Equations

DSolve[{eq1,eq2,...}, {y1[x],y2[x],...}, x]    	solve a list of differential equations        

DSolve[{eq1,eq2,...}, {y1,y2,...}, x]   		same task, pure functions as result

The following system:             		             
                                                                x'  -   y  =  t, 	 y'  - 4 x  = - 2
is solved:

Clear[x,y,f,g];
sy = {x'[t] - y[t] == t, y'[t] - 4 x[t] == - 2}

DSolve[sy, {x[t],y[t]}, t]//Flatten

Clear[x,y,f,g];
sy = {x'[t] - y[t] == t, y'[t] - 4 x[t] == - 2, 
              x[0] == 1,           y[0] ==  1   }

DSolve[sy, {x[t],y[t]}, t]

DSolve[sy, {x,y}, t]

11.1.4  Motion of Charged Particle in Static Crossed Homogeneous 
	    Electric and Magnetic Fields

D.M. Cook et al., Comp. in Phys., 6 (#5, Sep./Oct.92), 532f.

Equation of motion:

((m  (d^2)r)/(dt^2))  =  q (( dr)/dt) ×  B +    q  E

Initial data:

E   =   (0, E0, 0), B   =   (0, 0, B0) ;       r(0)  =  (0,0,0),     dr(0)/dt     =   (vx, vy, vz).

11.1.4.1   Setting up the Equations of Motion and Solving Them

Clear[m,r,v,a,B0, E0]
Efield = {0, E0, 0};   Bfield = {0, 0, B0};

r = { x[t], y[t], z[t]}; v = D[r,t]; a = D[v,t];

Force = q Efield + q Cross[ v, Bfield ]

eq = Table[ Force[[i]] == m a[[i]], {i,3}]

id = { x[0] == 0, x'[0] == vx, y[0] == 0, y'[0] == vy,
       z[0] == 0, z'[0] == vz };
       
Soln = DSolve[ Union[eq,id], {x,y,z}, t ]//Flatten;TableForm[Soln]

11.1.4.2   Checking the Solution to fulfil the euations of motion and the initial data

Simplify[eq /. Soln]

Simplify[{x[0], y[0], z[0]} /. Soln]

Simplify[{x'[0], y'[0], z'[0]} /. Soln]

11.1.4.3  Getting Expressions

X[t_] := Expand[ x[t] /. Soln] 
Y[t_] := Expand[ y[t] /. Soln] 
Z[t_] := Expand[ z[t] /. Soln]

rs = {X[t], Y[t], Z[t]} /. {m -> q B0/ω,  E0 -> vd B0};
TableForm[rs]

rf = rs /. {ω -> 1, vd -> 0.04, vx -> 0.05,
	vy -> .033, vz -> 0.03}; TableForm[rf]

11.1.4.3  Drawing Curves

pleb1 = ParametricPlot3D[ rf1, {t,0,20}, Boxed -> False,
 AxesLabel -> {"x","y","z"},
BoxRatios -> {1,1,1},ViewPoint->{-1.010, -2.400, 2.000}];

rf2 = rs /. {ω -> 1, vd -> 0.02, vx -> 0.05,
	vy -> .06, vz -> 0.03} ;
pleb2 = ParametricPlot3D[rf2, {t,0,20},Boxed -> False,
 AxesLabel -> {"x","y","z"}, 
BoxRatios -> {1,1,1},ViewPoint->{-1.010, -2.400, 2.000}];

Show[GraphicsRow[{pleb1,pleb2}]]

11.1.5	Changing the Dependent Variable of a Differential Equation

A differential equation can be transformed by changing the dependent variable   y(x)   of this equation. In most cases this is done by assuming that the new dependent variable is  z(x)  =  y(x)/w(x)  with the function  w(x)  chosen such that the new differential equation has a form more suitable for a given purpose. For example, in a second order differential equation one may remove the first derivative.  Or the new differential equation has standard functions as solutions.  Here, in a specific example, it is shown how this task may be accomplished. We start from Bessel's equation and transform it into a new equation without the first derivative.

Clear[m,w,x,z,y]
oldeq = y''[x] + y'[x]/x + y[x](1 - m^2/x^2)

Length[oldeq]

FullForm[oldeq]

cold = 
Table[Coefficient[oldeq, Derivative[k][y][x]], {k,0,2}]

The above list gives the coefficients of the various derivatives. This is used  to find a new differential equation, which is the old differential equations with the product  y(x) = z(x) w(x)  inserted for the old dependent variable. The resulting equation is rearranged to give better insight.

neweq = 
Sum[ D[z[x] w[x], {x,k}] cold[[k+1]], {k,0,2}]//Expand

cnew = 
Table[Coefficient[neweq, Derivative[k][z][x]], {k,0,2}]

fneweq = Sum[D[z[x], {x,k}] cnew[[k+1]], {k,0,2}]

Now  we prescribe the new equation to have no first derivative. This yields a first order differential equation for the unknown function  w(x);  this is solved and the resulting solution is inserted into the new differential equation. For this purpose it is advantageous to get the solution for w(x) in a pure function form
(cf.  11.1.2   and  21.3 ).

heq = cnew[[2]] == 0

s1 = DSolve[heq, w, x]//Flatten

fnew = fneweq /. s1

normf = Coefficient[fnew, Derivative[2][z][x]]

fnew = fnew / normf //Expand

fdeq = Collect[fnew, z[x]/x^2] == 0

11.2   Numerical Solution of Differential Equations

In Mathematica it is quite easy to get z solutions of ordinary differential equations; and to plot the solutions.

More information on the corresponding package may be found at:
http://reference.wolfram.com/mathematica/tutorial/NDSolvePackages.html

NDSolve[{eq1, eq2, ...}, y, {x, xmin, xmax}]                                                            
	 find a numerical solution for the function  y  with x  in the range  xmin to xmax      
                                                           
NDSolve[{eq1, eq2, ...}, {y1, y2,...}, {x, xmin, xmax}]                                               
               find numerical solutions for several functions   yi

In the following is assumed that the numeric solution is labelled  solution.

sol[x_] = y[x] /. solution    use the list of rules for the function y to get a function sol[x] 					  for the solution y[x]

Plot[y[x] /. solution, {x,xmin, xmax} ]    plot the solution to a differential equation

11.2.1  Numeric Solutions of Systems of Differential Equations 
	for One Dependent Variable

11.2.1.1  Example 1

Clear[y,x]; xmin = 0.; xmax = 6;
s = NDSolve[ {y''[x] == - y[x], y[0] == 0, y'[0] == 1}, 
			  							   y,{x,xmin,xmax}]//Flatten

y[2.31]    /. s
y'[2.31]   /. s
y''[2.31]  /. s

np = 5; xk := k (xmax - xmin)/np; 
ts = Table[{k,N[xk],y[xk] /. s,y'[xk] /. s },  {k,0,np}]//Chop;

TableForm["k    x     y(x)         y'(x)  "]
TableForm[ts]

py=Plot[y[x]/.s,{x,xmin,xmax},AxesLabel->{x,y[x]}];

dpy=Plot[y'[x]/.s,{x,xmin,xmax},AxesLabel->{x,y'[x]}];

Show[GraphicsRow[{py,dpy}]]

11.2.1.2  The function defined by the numeric solution of the differential equation

It is possible to use the numeric solution of a differential equation for defining a function, which may be used for getting values of the solution or plots. As an example we take a nonlinear differential equation with initial  data:

Clear[sol];
sol[t_]=First[x[t]/.NDSolve[{x''[t]==(Sqrt[x[t]]) x'[t] + 1,x[1]==7,x'[1]==-1},{x[t]},{t,1,2},MaxSteps->10^4]]

sol[1.5]

Plot[sol[t],{t,1,2}, ImageSize -> 250]

11.2.1.3  Motion of mass point in gravity field with and without Newtonian friction.

Clear[x,y,z,r,v,b]
r[t_] = {x[t], y[t]};
v[t_] = D[r[t],t];
b[t_] = D[v[t],t];
m = 1; g = 10; a = .3; 
sys = Thread[ m b[t] == {0, -m g}   ]

sysa = m b[t] == 
{0, -m g} - a v[t] Sqrt[ v[t] . v[t] ] //Thread

anf = {x[0] == 0, y[0] == 0, x'[0] == 2, y'[0] ==10};

sol = NDSolve[ Join[sys,anf], {x,y}, {t,0,4}] //Flatten ;

sola = NDSolve[ Join[sysa,anf], {x,y},{t,0,2}] //Flatten ;

p  =ParametricPlot[{x[t],y[t]}/.sol,{t,0,4}];
pa=ParametricPlot[{x[t],y[t]}/.sola,{t,0,2},  
         PlotStyle->Dashing[{0.01`}]];

Show[p,pa,AxesLabel->{"x","y"}, AspectRatio ->   0.4]

We assume the plane y = 0 to be the ground. So the trajectories should end there. We could find the appropriate value of the independent variable  t   by trial and  error. A more systematic approach uses the command  FindRoot[]  to find these values of   t.  Note that one must use this command with two starting values.

te = t/.FindRoot[y[t] /. sol   == 0., {t,2,4}]
ta = t/.FindRoot[y[t] /. sola  == 0., {t,1,2}]

p  =ParametricPlot[{x[t],y[t]}/.sol,{t,0,te}];
pa=ParametricPlot[{x[t],y[t]}/.sola,{t,0,ta}, 
  PlotStyle->Dashing[{0.01`}]];

Show[p,pa,AxesLabel->{"x","y"}, AspectRatio -> Automatic, ImageSize -> 200]

11.2.2  Numeric Solutions of Systems of Differential Equations
	for Several Dependent Variables.

sy = { -(y'[t]) + x''[t] == 0,  x'[t] + y''[t] 
   == 0,  z''[t] == 0,  x[0] == 0, y[0] == 0,
   z[0] == 0, x'[0] == 2, y'[0] == 0, z'[0] == 1/3 }

tmax = 11;
sosy = NDSolve[sy, {x,y,z}, {t, 0, tmax}]

ParametricPlot3D[{x[t],y[t],z[t]}/.sosy,{t,0,tmax},BoxRatios->{1,1,1.5`},AxesLabel->{"x","y","z"}, ImageSize -> 200]

11.2.3   Options for NDSolve

??NDSolve

sy = { -(y'[t]) + x''[t] == 0,  x'[t] + y''[t] 
   == 0,  z''[t] == 0,  x[0] == 0, y[0] == 0,
   z[0] == 0, x'[0] == 1, y'[0] == 0, z'[0] == 1/3 };

NDSolve[sy, {x,y,z}, {t,0,100}, MaxSteps -> 200]

NDSolve::mxst :  "Maximum number of 200  steps reached at the point t ==  28.340042682535987`. >>

So one must increase the number of steps. Thereafter everything works properly:
As stressed before: in general, it is not necessary to choose the MaxSteps number.

NDSolve[ sy, {x,y,z}, {t,0,100}, MaxSteps -> 1000]

??AccuracyGoal

??PrecisionGoal

??WorkingPrecision

??StartingStepSize

The last four options may be very important for achieving sufficient accuracy of the solutions obtained by numeric  integration of differntial equations. This is  discussed in the next  section.

??Method

Various methods of solution can be selected. Find details via Help/Find Selected/Functions/NDSolve .
See also 11.2.8

11.2.4  Checking  the Accuracy of Numeric Solutions of Differential Equations.

A method to check the accuracy of numeric solutions of differential equations  is to run  NDSolve[]  at least twice with different values assigned to the options AccuracyGoal, PrecisionGoal  and  WorkingPrecision (s. 11.2.3). Another method for such a check applicable to differential equations describing a conservative system is to compare the initial and final values of total energy. The latter check, however, is a rather weak one. A third check is to run the 
system in a second run in the reverse direction and to check whether the final data of the latter run agree with the initial date of the first run. All this is shown for the nonlinear system discussed  in more detail in 11.2.6 .

Equations of Motion

Clear[t,r,x,y,vx,vy,sys,en,tmax, force]
r[t_] = {x[t], y[t]};  v[t_] = {vx[t],vy[t]};
force = {-x[t] - 2 x[t] y[t], -y[t] + y[t]^2 - x[t]^2} ;
sys = D[Join[r[t],v[t]],t]  == Join[v[t], force ]//Thread

en = (vx[t]^2 + vy[t]^2)/2 + (x[t]^2 + y[t]^2)/2 + 
	  x[t]^2 y[t] - y[t]^3/3;

ad = {x[0], y[0], vx[0], vy[0]};

Run with default precision for a short time

anf =  {0, 0.1, 0.15, 0.55}; tmax = 15.;

sol = NDSolve[ Join[sys,Thread[ad== anf]], {x,y,vx,vy}, {t,0,tmax}]//Flatten ;

eni = ( en /. sol ) /. t -> 0

enf =( en /. sol ) /. t -> tmax

(enf - eni)/eni

This check is useful inasmuch as a large difference in the values of  the initial and the final energies indicates a serious inaccuracy. But energy is a stationary quantity; so it is insensitive to errors in the data used to evaluate it.  Therefore this check represents more a necessary than a sufficient condition.

A more sensitive check is to evaluate the coordinates and velocities at the final time point,  to run the system backwards in time till the initial time and to compare the final coordinates and velocities of this run to the initial coordinates and velocities used in the first run. The time reversal just described requires to use the final coordinates and the negative values of the final velocities of the first run as initial data for the second run.

Return with default precision for the same short time

{x[tmax],y[tmax],vx[tmax],vy[tmax]} /. sol//Flatten

anfr = % {1,1,-1,-1}

sor =  NDSolve[ Join[sys,Thread[ad== anfr]], {x,y,vx,vy}, {t,0,tmax}]//Flatten ;

{x[tmax],y[tmax],vx[tmax],vy[tmax]} /. sor//Flatten

anf

Now the final data of the reversed run  can be compared to the inital data  (anf). Agreement is quite good. The different sign of the velocity is correct as in the second run the trajectory returns to the starting point of the first run. The trajectories  in coordinate space are indistinguishable.

p1=ParametricPlot[r[t]/.sol,{t,0,tmax},AxesLabel->{"x","y"},PlotLabel->"First run", AspectRatio ->  1];
p2=ParametricPlot[r[t]/.sor,{t,0,tmax},AxesLabel->{"x","y"},PlotLabel->"Return run", AspectRatio ->  1];
Show[GraphicsRow[{p1,p2}]]

The solutions of this highly nonlinear system near the stability limit require appropriate options for runs over long time intervalls; without these the solutions become inaccurate soon.  But these options increase the CP time, the running time and data storage space very much.  The accuracy is checked  by running backwards from the end point to compare these end data with the initial data.

Run with default precision for a long time

tmax = 150.;

{time,sol0} = NDSolve[ Join[sys,Thread[ad== anf]], {x,y,vx,vy}, {t,0,tmax},
MaxSteps -> 5000]//Flatten//Timing

g0x = Plot[ x[t]/.sol0, {t,0,tmax}, 
		   AxesLabel -> {"t","x[t]"},
           PlotStyle -> {{AbsoluteThickness[1],Dashing[{.01}]}},
           PlotPoints -> 500];
g0y = Plot[ y[t]/.sol0 , {t,0,tmax}, 
		   AxesLabel -> {"t","y[t]"},
           PlotStyle -> {{AbsoluteThickness[1],Dashing[{.01}]}},
           PlotPoints -> 500];

Return with default precision for a long time

{ x[tmax], y[tmax], vx[tmax], vy[tmax] } /. sol0

anfr = % {1,1,-1,-1}

{time,sor0} = NDSolve[Join[sys,Thread[ad==anfr]], {x,y,vx,vy}, {t,0,tmax}, 
	                  MaxSteps->5000]//Flatten//Timing

{ x[tmax], y[tmax], vx[tmax], vy[tmax] } /. sor0

anf

eni = (en /. sol0)  /.  t -> 0

enr = (en /. sor0)  /.  t -> tmax

It is obvious that the end values of the return run differ strongly from the intial data  (anf) of the first run, although the various values of total energy do not differ very much.  To get reliable values for the solutions over the full time intervall appropriate options must be prescribed in NDSolve[].

Run with high precision for a long time

{time,sol1} = NDSolve[Join[sys,Thread[ad==anf]], {x,y,vx,vy}, {t,0,tmax}, 
	          AccuracyGoal  -> 20,  WorkingPrecision -> 25,
	          PrecisionGoal -> 20, MaxSteps -> 400000]//Flatten//Timing

The above warning tells us that the statements describing the differential equations and the initial conditions are not given with the working precision (25).  The final tests show that the precision obtained suffices for our purposes.

g1x = Plot[x[t]/.sol1, {t,0,tmax}, 
		   AxesLabel -> {"t","x[t]"},
           PlotStyle -> AbsoluteThickness[1],
           PlotPoints -> 500];
g1y = Plot[ y[t]/.sol1, {t,0,tmax}, 
		   AxesLabel -> {"t","y[t]"},
           PlotStyle -> AbsoluteThickness[1],
           PlotPoints -> 500];

Return with high precision for a long time

SetPrecision[{ x[tmax], y[tmax], vx[tmax], vy[tmax] } /. sol1,30]

anfr = % {1,1,-1,-1}

{time,sor1} = NDSolve[Join[sys,Thread[ad==anfr]], {x,y,vx,vy}, {t,0,tmax}, 
	          AccuracyGoal  -> 20, WorkingPrecision -> 25,
	          PrecisionGoal -> 20, MaxSteps->400000]//Flatten//Timing

{ x[tmax], y[tmax], vx[tmax], vy[tmax] } /. sor1

anf

eni = ( en /. sol1 )  /. t -> 0

enr = ( en /. sor1 )  /. t -> tmax

Now the end values of the return run agree well with the initial data of the first run. So one can be sure that the solutions are accurate over the whole intervall of the integration. Below the coordinates of the points are compared as obtained with the stronger options (continuous line) and with unsufficient options (dashed line).

g1x = Plot[x[t]/.sol1, {t,0,tmax}, 
		   AxesLabel -> {"t","x[t]"},
           PlotStyle -> {AbsoluteThickness[0.5],Black},
           PlotPoints -> 500];
g1y = Plot[y[t]/.sol1, {t,0,tmax}, 
		   AxesLabel -> {"t","y[t]"},
           PlotStyle -> AbsoluteThickness[0.5],
           PlotPoints -> 500];

Show[g1x,g0x, ImageSize -> 600,
PlotLabel -> "__  high precision,   --- default precision"]

Show[g0y,g1y,ImageSize->600,PlotLabel->"__  high precision,   --- default precision"]

Diagrams in x,y-space are less suitable to show the difference between accurate and inaccurate solutions. From the outputs above it is seen that the options for higher precision entail about a factor 100 in CP time. The CP times, however,  do not reflect the real running times, wich last still longer.  So other programming  languages compare favourably for such tasks.
In fact, a few years ago, using integration routines *) written in C++ and MATLAB graphics  software,  N. Ivanovic  needed about 2 to 3 sec computer time to prepare figures  corresponding to those given above. At that time Mathematica needed about 20 minutes to get the results.

*) W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery: Numerical Recipies in C.  2nd ed., 1992, Cambridge University Press.

11.2.4.1     Stiff  Differential  Equations

Now the difficulties occuring in the numeric solution of a stiff differential equation are discussed at an example.
These are notoriously difficult in numeric solutions.

sys = y'[t] == y[t] - t^2;    DSolve[sys,y[t],t]

This is the general solution. Note that it consists of a polynomial and of an exponential function. The latter grows much faster with increasing time as compared to the former; in general, it will win over the polynomial in the long run. But there are special initial data, which render the constant C[1] = 0. But in a numeric solution round-off errors mix it into the solution so that its initial suppression gets lost gradually.  This explains why both numeric solutions deviate from the corresponding analytic solution (displayed by dashed curves ).

anf1= y[0] == 2; so1 = DSolve[{sys,anf1},y[t],t]

anf2= y[0] == 2.00001 ; so2 = DSolve[{sys,anf2},y[t],t]

tmax=17.5`;
p1=Plot[y[t]/.so1,{t,0,tmax},PlotStyle->{{RGBColor[0.9`,0,0],Dashing[{0.02`}]}}] ;
p2=Plot[y[t]/.so2,{t,0,tmax},PlotStyle->{{RGBColor[0,0.7`,0],Dashing[{0.05`}]}}] ;
Show[p1,p2,PlotRange->All]
no1=Flatten[NDSolve[{sys,anf1},y,{t,0,tmax}]];
no2=Flatten[NDSolve[{sys,anf2},y,{t,0,tmax}]];
pn1=Plot[y[t]/.no1,{t,0,tmax},PlotStyle->RGBColor[0.9`,0,0]];
pn2=Plot[y[t]/.no2,{t,0,tmax},PlotStyle->RGBColor[0,0.7`,0]];
Show[pn1,p1,PlotRange->All]
Show[pn2,p2,PlotRange->All]

Show[p1,pn1,pn2,p2,PlotRange->All,AxesLabel->{"t","y(t)"},BaseStyle->{FontSize->10},PlotLabel->Row[{"Red = sol1,  Green = sol2,  Dashed = numeric\n\n"}]]

The routines Mathematica uses now are so good that the differences visible in earlier versions are no longer present.

11.2.5 	Phase Space Diagrams

The phase space of a dynamical system, whose degree of freedom is   f,  is a 2f-dimensional space. It is made up of the   f  pairs of dependent variables describing the state of the system. To each degree  there correponds a pair of variables, the position coordinate and the corresponding (generalized) velocity or the corresponding (canonical) momentum. The solutions obtained for x(t), x'(t), y(t), y'(t), ...give a curve in this 2f-dimensional space, the phase curve. A point on this curve, the phase point, describes the state of the system corresponding to the time  t. In this and the next section it is shown how one may draw phase space diagrams from numeric solutions of the equations of motion. In this section such diagrams for  f = 1 are treated.  If  f > 1  then the dimension of phase space is too large to permit one to draw the full space. One is restricted toprojections. An important representation of this type is the Poincaré Map  treated in the next section.

11.2.5.1  Damped Harmonic Oscillator

deq = x''[t] + .3 x'[t] + 5 x[t] ; tmin = 0; tmax = 20;
dosc = NDSolve[ {deq == 0, x[0]==0, x'[0] == 1}, x, {t,tmin,tmax}]

ParametricPlot[{x[t],x'[t]}/.dosc,{t,tmin,tmax},ImageSize -> 150,PlotRange->{{-0.45`,0.45`},{-1,1}},AxesLabel->{x,x'},Ticks->None]

11.2.5.2  Damped Mathematical Pendulum

deq = p''[t] +  p'[t] + 5 Sin[p[t]] ; tmin = 0; tmax = 6;
dap = NDSolve[ {deq == 0, p[0]==0, p'[0] == 8}, p, {t,tmin,tmax}];

dapp=Plot[p[t]/.dap,{t,tmin,tmax},AxesLabel->{t,"ϕ(t)"},PlotRange->{0,8},Ticks->{{1,2,3,4,5,6},{0,π,2 π}},PlotLabel->"Damped Pendulum"];
dappp=Plot[p'[t]/.dap,{t,tmin,tmax},AxesLabel->{t,"ϕ'(t)"},PlotRange->{-2,8},Ticks->{{1,2,3,4,5,6},{1,2,3,4,5,6,7,8}},PlotLabel->"Damped Pendulum"];

Show[GraphicsRow[{dapp,dappp}], ImageSize ->   550]

plc=ParametricPlot[{p[t],p'[t]}/.dap,{t,tmin,tmax},AxesLabel->{"ϕ","ϕ'"},Ticks->{{0,π,2 π},{-2,0,2,4,6,8}},PlotStyle->{GrayLevel[0],Thickness[0.004]},PlotRange->{{0,8},{-2,8}}, ImageSize -> 180 ]

The first two figures show the angle and the angular velocity as functions of time.  The full phase curve is shown just before.  

But the phase space diagram does not contain information on time. This can be done by adding points on the phase curve corresponding to equal time steps. The method to achieve this is shown at first for a small number of intervalls.

np=5;tk:=((k (tmax-tmin))/np);Print[{"k"," time","      ϕ  ","   ϕ'      "}] Do[Print[Column[{{k,N[tk+0.001`],{p[tk],p'[tk]}/.dap}},Center]],{k,0,np}]

np = 20; tk := k (tmax - tmin)/np; 
lp = Table[Flatten[{p[tk],p'[tk]} /. dap], {k,0,np}];

plp=ListPlot[lp,Ticks->{{0,π,2 π},{-2,0,2,4,6,8}}];

plt=Show[plc,plp,Prolog->PointSize[0.01`],PlotRegion->{{0.01`,0.99`},{0.01`,0.99`}}];Show[GraphicsRow[{plp,plt}],ImageSize->550]

11.2.6	Poincaré  Maps

The phase space of a system of two degrees of freedom is 4-dimensional; so there is no means to represent it in full. The situation is even worse if f > 2. One must recur to projections of the phase space. A particularly important representation is the Poincaré Map. The phase curve (cf. preceeding section) intersects an appropriately chosen plane (a subspace of the full phase space) frequently. Each intersection is marked by a dot. This pattern of dots is called a Poincaré Map. It displays important information on a dynamical system. If there are isolating integrals of the motion (as e.g. energy or angular momentum) these will restrict the subspace available to the phase curve considerrably. Some of such integrals may not be valid globally (.i.e in the whole phase space) but only locally (in a restricted subdomain of phase space). The action of such local integrals may be rendered visible indirectly in Poincaré Maps. They lead to rather regular shapes of the trajectories in coordinate and in phase space and to regular patterns in the Poincaré Maps; whereas the lack of such integrals leads to complex phase curves and to chaotic distributions of the points in the map.

This is shown with a specific example developed by the French astronomer Hénon describing the motion of a star of unit mass in the field of the other stars making up a galaxis. The equations of motion are 2-dimensional; the force on the star is nonlinear		  		
				F   =   ( - x  - 2 x y,  - y  +  y2  -  x 2 ).
There is a potential  V(x,y); so  energy is conserved. 

	V{x,y}   =   (x2  +  y2 )/2   +  x2  y  -  y3 /3 ;      E   =   ((OverDot[x])2  + (OverDot[y]) 2)/2  +  V(x,y)  =  const.

For small values of total energy there exists an additional local integral of the motion. At first the equations of motion are solved by NDSolve. It is advantageous to represent the equations of motions as a system of first order differential equations.  For, on some machines the solution process runs faster than that for the system of second order equations.  In any case the velocities are needed for the map. Application of  NDSolve gives an (internal) representation of the functions   x(t), y(t), vx(t), vy(t).  This can be used to calculate the tranjectory in coordinate space or phase space. For 
the Poincaré map the instants must be found for which the phase curve intersects the plane (say x = 0) selected as the subspace of the phase space for the particular map. Each of these times inserted into the solution gives a dot in the  y, vy-plane.

Experience (s. end of 11.2.4.) shows that Mathemtica is not very well suited for preparing Poincaré maps for non-linear systems displaying chaotic behavior.  It requires too much time and space on central memory and disks for calculations giving reliable results. So the examples given below should be regarded as crude illustrations.

The trajectorys in x,y-space, x(t) and y(t), finally the Poincaré map, will be shown for two cases:   
	1. For a small energy eno =  .01, where motion is still well ordered;   
	2. for an energy enc  =  .161, which is slightly below the limit  enl = 1/6 = .166666...  for 			bounded motion,  where motion is chaotic.  
The potential of this system is shown in perspective and by a contour drawing below.  V = 1/6 corresponds to the three saddle points.

Plotting the potential

v[x_, y_] = (x^2 + y^2)/2 + x^2 y - y^3 / 3;

plv=Plot3D[v[x,y],{x,-1.2,1.2`},{y,-1`,1.05`},PlotPoints->30,PlotRange->{0,0.2`},AxesLabel->{"x","y","   V(x,z)"},ViewPoint->{1,-0.6`,1},Ticks->{Automatic,Automatic,{0.`,0.1`,0.2`}}];
plc=ContourPlot[v[x,y],{x,-1.2,1.2`},{y,-1.2`,1.05`},ColorFunction->Hue];

Show[GraphicsRow[{plv,plc}],ImageSize->500]

Equations of motion in Hamiltonian form

r[t_] = {x[t], y[t]}; v[t_] = {vx[t],vy[t]};
force =
	{-x[t] - 2 x[t] y[t], -y[t] + y[t]^2 - x[t]^2} ;		  
sys =D[Join[r[t],v[t]], t] == Join[v[t],force] //Thread

en = (vx[t]^2 + vy[t]^2)/2 + (x[t]^2 + y[t]^2)/2 + 
	  x[t]^2 y[t] - y[t]^3/3;

Orderly motion

At first the initial data are chosen such that the mortion is orderly.

x0 = 0.;  y0 = .01;  vx0 = .141;  vy0 = 0.; 
tmax = 251 Pi//N

eno = en /. 
  {x[t] -> x0, y[t] -> y0, vx[t] -> vx0, vy[t] -> vy0}

anf={x[0]==x0,y[0]==y0,vx[0]==vx0,vy[0]==vy0}

solo=NDSolve[Join[sys,anf],{x,y,vx,vy},{t,0,tmax},MaxSteps->6500]//Flatten

solo2=NDSolve[Join[sys,anf],{x,y,vx,vy},{t,0,2tmax},MaxSteps->15000]//Flatten

The nonlinear terms lead to a nonclosed trajectory which rotates slowly.The ellipse of linear motion is shifted and slightly deformed. The figure below showing the trajectory for a longer time intervall displays  the threefold symmetry of the potential well.

ParametricPlot[r[t]/.solo,{t,0,6.5`},AxesLabel->{"x","y"},AspectRatio->Automatic,PlotStyle->{AbsoluteThickness[1`],GrayLevel[0]}, ImageSize ->  500]

ParametricPlot[r[t]/.solo,{t,0,75.`},AxesLabel->{"x","y"},AspectRatio->Automatic,PlotStyle->{AbsoluteThickness[1`],GrayLevel[0]}, ImageSize ->  500]

ParametricPlot[r[t]/.solo,{t,0,tmax},AxesLabel->{"x","y"},AspectRatio->Automatic,PlotStyle->AbsoluteThickness[1],PlotPoints->3000, ImageSize ->  500]

The coordinates as functions of time show also a rather regular behavior. They show oscillations with a regular modulation of the amplitude.

Plot[x[t]/.solo,{t,0,150.`},PlotPoints->2500,PlotStyle->AbsoluteThickness[1],AxesLabel->{"t","x"}]

Plot[y[t]/.solo,{t,0,150.`},PlotPoints->2500,PlotStyle->AbsoluteThickness[1],AxesLabel->{"t","y"}]

Plot[x[t]/.solo,{t,0,tmax},PlotPoints->3500,PlotStyle->{Thickness[0.001`],GrayLevel[0]}, ImageSize ->  500,AxesLabel->{"t","x"}]

Plot[y[t]/.solo],{t,0,tmax},PlotPoints->3500,PlotStyle->{Thickness[0.001`],GrayLevel[0]}, ImageSize ->  500, AxesLabel->{"t","y"}]

Chaotic motion

Now the trajectory of chaotic motion is computed. In order to get reliable data the various precision options must be increased as discussed in 11.2.4.  The trajectory is rather irregular.  The same applies to the coordinates as function of time as shown in the last two figures of 11.2.4.

x0=0.`; y0=.1`; vx0=.15`; vy0=0.54`;
anf={x[0]==x0,y[0]==y0,vx[0]==vx0,vy[0]==vy0}

enc=en/.{x[t]->x0,y[t]->y0,vx[t]->vx0,vy[t]->vy0}

solc=NDSolve[Join[sys,anf],{x,y,vx, vy},{t,0,tmax},AccuracyGoal->15,WorkingPrecision->25,PrecisionGoal->15,MaxSteps->26250]//Flatten

solc2=NDSolve[Join[sys,anf],{x,y,vx, vy},{t,0,2tmax},AccuracyGoal->15,WorkingPrecision->25,PrecisionGoal->15,MaxSteps->50250]//Flatten

ParametricPlot[{x[t],y[t]}/.solc,{t,0,150.`},AxesLabel->{x,y},PlotStyle->AbsoluteThickness[1],AspectRatio->Automatic]

ParametricPlot[{x[t],y[t]}/.solc,{t,0,tmax},AxesLabel->{x,y},PlotStyle->AbsoluteThickness[0.5],PlotPoints->15000,AspectRatio->Automatic]

For drawing the Poincar› map the times  ti  must be found for which   x(ti) = 0. This is done with the help of  FindRoot. In order to get two starting values for this command, at first a list of values  x(tk)   is computed. The difference  tk+1 -   tk   must be chosen small enough; it can be read off from the figures showing  x(t)  given above.  This list is scanned for changes of sign;  two adjacent values having different signs are used as starting values for FindRoot.This gives a list of the  ti .Then a list of values  { y(ti), vy(ti) }  is computed and plotted. At first ths is done for a short running time to show the working of the method.

solo2[[1]]

fx[t_]= x[t]/.solo2[[1]];
fy[t_]= y[t]/.solo2[[2]];
fpx[t_]= vx[t]/.solo2[[3]];
fpy[t_]= vy[t]/.solo2[[4]];

Presenting the method of finding roots for a short running time

pas = 1;
tx = Table[fx[t], {t,0,10,pas}]//Chop

The changes of sign are found be multiplying adjacent values and evaluating the sign of this product:

erstes und letztes element der liste wird weggeschnitten, anschl. elementweise multiplikation, signum-Fkt, negativer wert bei vorzeichenwechsel

Sign[Drop[Chop[tx],1] Drop[Chop[tx],-1]]

The positions of -1 in the list above give the positions of the value of x preceeding a change in sign.
This is used as a starting value for FindRoot[] searching for the value  (Subscript[t, i]) of the zero x[(Subscript[t, i])] = 0.

li=Flatten[Position[
Sign[Drop[tx,1] Drop[tx,-1]],-1]]

Table[FindRoot[fx[t] == 0.,{t,li[[k]]}],
{k,Length[li]}]

An approximate value t0 of the zero is given by the point where the straight line joining the
neighbouring points  ((Subscript[t, i]), (Subscript[x, i])) and  ((Subscript[t, i+1]), (Subscript[x, i+1])) crosses the t-axis. The corresponding formula
is shown with  (Subscript[t, i])= t1,  (Subscript[x, i])= x1 and  (Subscript[t, i+1])= t2, (Subscript[x, i+1])= x2 :

t0 = t1 -  x1 ((t2 - t1)/(x2 - x1)) =  ((t1 x2 - x1 t2)/(x2 - x1))

tax =
(pas*(li-1)*fx[pas*li]-pas*li*fx[pas*(li-1)])/
(fx[pas*li]-fx[pas*(li-1)]) ;

fx[tax]

Due to the large step width  pas  the approximate values of the zeros (Subscript[t, i])given in the list  tax
are still somwhat crude and the corresponding values x((Subscript[t, i])) are small, but not very small.
But this improves with a smaller step width.

pas = .1;
tx = Table[fx[t], {t,0,10,pas}]//Chop ;

li=Flatten[Position[
Sign[Drop[tx,1] Drop[tx,-1]],-1]];

tax =
(pas*(li-1)*fx[pas*li]-pas*li*fx[pas*(li-1)])/
(fx[pas*li]-fx[pas*(li-1)])

fx[tax]

{fy[tax],fpy[tax]}//Transpose

The list above gives the pairs (y, vy) for the three times given in tax.

Another method for finding the zeros  (Subscript[t, i]) has been given by C. Berger:

tx = Table[x[t] /. solo /. t -> k, {k,0,10}];
ts = Module[{l = {}},  Do[ If[ tx[[k]] tx[[k+1]] <= 0, l = 
	   		Append[l, (FindRoot[ x[t] /. solo == 0., {t,k,k-1}])] ],{k,Length[tx]-1} ]; l]

Drawing the Poincaré map in the  y, py-plane for oderly motion

pas = .1;
tx = Table[fx[t], {t,0,2tmax,pas}]//Chop ;

li=Flatten[Position[
Sign[Drop[tx,1] Drop[tx,-1]],-1]];

tax =
(pas*(li-1)*fx[pas*li]-pas*li*fx[pas*(li-1)])/
(fx[pas*li]-fx[pas*(li-1)]) ;

Max[Abs[fx[tax]]]

py=ListPlot[Transpose[{fy[tax],fpy[tax]}],Axes->None,RotateLabel->False,Frame->True,FrameLabel->{"y","vy"}]

Drawing the Poincaré map in the  y, py-plane for chaotic motion

fx[t_]= x[t]/.solc2[[1]];
fy[t_]= y[t]/.solc2[[2]];
fpx[t_]= vx[t]/.solc2[[3]];
fpy[t_]= vy[t]/.solc2[[4]];

pas = .1;
tx = Table[fx[t], {t,0,2tmax,pas}]//Chop ;

li=Flatten[Position[
Sign[Drop[tx,1] Drop[tx,-1]],-1]];

tax =
(pas*(li-1)*fx[pas*li]-pas*li*fx[pas*(li-1)])/
(fx[pas*li]-fx[pas*(li-1)]) ;

Max[Abs[fx[tax]]]//ScientificForm

py=ListPlot[Transpose[{fy[tax],fpy[tax]}],Axes->None,RotateLabel->False,Frame->True,FrameLabel->{"y","vy"}]

11.2.7 	Finding Eigenvalues by the Shooting Method

The shooting method is a means to compute eigenvalues of a boundary problem defined by a one-dimensional differential equation and boundary conditions fixed at the end of the interval whereinthe equation is defined. 
At first a trial value is inserted for the eigenvalue and the differential equationis solved by numeric integration for an initial problem with initial data fulfilling the boundary condition(s) at one end. The integration goes as far as the other end. In most cases inspection of the solution shows that there,the boundary condition(s) is (are) not fulfilled. But the result gives a hint how to improve the trial value for the eigenvalue. This procedure is repeated till a solution fulfilling all boundary conditions is obtained. The trial value then used is an approximation to an eigenvalue of the boundary value problem. 

This method is displayed for the eigenvalue problem of the one-dimensional harmonic oscillator in quantum mechanics. The differential equation expressed in dimensionless variables is :

(((∂^2)y)/(∂(x^2)))    +   (2 ν+1  -(x^2)) y  =  0

ν  is the eigenvalue parameter; it is related to energy  E  by   E =  ℏω (2 ν + 1)/2.  The boundary condtions are:

x = ± ∞:  y = 0.

Analytic solution leads to   E/ℏω  =  1/2, 3/2, ...,  i.e. ν = 0, 1, ... ; see § 11.1.1 This is now shown by the shooting method. In place of  x = ± ∞  a value of sufficient magnitude is chosen; experience shows that for low  ν  x = 5  suffices; the value of the first derivative must be small, the specific value is uncritical. The end value xe is determined empirically such that the resulting curve fits into the picture.

Clear[x,y,sh,nu,xe]

11.2.7.1	  Showing the shooting method in detail for  n = 0

sh[nu_,xe_] :=
NDSolve[ {y''[x]  + (2 nu + 1 - x^2) y[x] == 0, 
	y[-5] == 4 10^-6, y'[-5] == 5 10^-5}, y[x], {x,-5,xe}]

SetOptions[Plot, PlotStyle -> Thickness[.005]];
dd = PlotStyle -> Dashing[{.01}];
dt = PlotStyle -> Dashing[{.02,.01,.0025,.01}];

so= Table[0,{5}]; pp =so;

so[[1]]=sh[-1.5`,xe=-0.5`];pp[[1]]=Plot[y[x]/.so[[1]],{x,-5,xe},Evaluate[dd]];so[[2]]=sh[-0.5`,xe=1];pp[[2]]=Plot[y[x]/.so[[2]],{x,-5,xe},Evaluate[dd]];so[[3]]=sh[-0.01`,xe=6];pp[[3]]=Plot[y[x]/.so[[3]],{x,-5,xe},Evaluate[dd]];so[[4]]=sh[0,xe=7];pp[[4]]=Plot[y[x]/.so[[4]],{x,-5,xe}];so[[5]]=sh[0.01`,xe=6];pp[[5]]=Plot[y[x]/.so[[5]],{x,-5,xe},Evaluate[dt]];

Show[pp[[1]],pp[[2]],pp[[3]],pp[[4]],pp[[5]],PlotRange->{-5,5},AxesLabel->{"x","y(x)"},PlotLabel->"QM Harmomiic Oscillator:  ν =  -1.5, -.5, -.1, 0, .1 . \n\n",PlotRegion->{{0.01`,0.99`},{0.01`,0.99`}},
Epilog -> {Text["ν = -1.5",{-2.1,4.75}],Text[-0.5,{0.19,4.75}],
Text["ν = - .1",{2.47,4.75}],Text["ν = 0.1",{2.47,-4.75}],
Text["ν = 0 ",{5,0.55}]}]

The boundary conditions are best approximated by the continuous curve corresponding to  ν = 0. Note the different behaviour  of the solutions belonging to  ν < 0 (dashed curves) or ν > 0 (dashed-dotted curves). 
Since the curve for   ν = 0  has no node in the interior of the interval, the eigenvalue must be the lowest one.

11.2.7.2  Automating the routine for finding the eigenvalues

Now the method is automated so that a routine searches for the characteristic values. A function is defined to get the end value y(xe)  as a function of  ν =nu. This function is used in FindRoot[]. In order to avoid troubles with FindRoot[] one must either 
	1) apply the operator Hold[] to the definition of tr[n_u]. This holds the evaluatin till te[nu] is 	     used in FindRoot[]. 
	     This hold is released by the command Release.  or
	2) restrict the arguments of the function  te[] to numeric values as done at the end of this 		    subsubsection.

xe = 5;

tet[nu_] := Hold[y[xe] /. NDSolve[ {y''[x]  + (2 nu + 1 - x^2) y[x] == 0, 
	y[-5] == 4 10^-6, y'[-5] == 5 10^-5}, y, {x,-5,xe}]]

FindRoot[tet[nnn],{nnn,.8,1.3}]
tet[nnn /. %]//Release

FindRoot[tet[nnn],{nnn,1.8,2.3}]
tet[nnn /. %]//Release

FindRoot[tet[nnn],{nnn,2.8,3.3}]
tet[nnn /. %]//Release

11.2.7.3   The variant wherein  nu  is restricted to numeric values

ten[nu_?NumericQ] := y[xe] /. NDSolve[ {y''[x]  + (2 nu + 1 - x^2) y[x] == 0, 
	y[-5] == 4 10^-6, y'[-5] == 5 10^-5}, y, {x,-5,xe}]

FindRoot[ten[nnn],{nnn,.8,1.3}]
ten[nnn /. %]

FindRoot[ten[nnn],{nnn,1.8,2.3}]
ten[nnn /. %]

11.3 Solutions of Partial Differential Equations.

11.3.1  Analytical Solution

Remember that while an ordinary  differential equation determines the solution up to some integration constants the solution of a partial differential equation contains arbitray functions. These are denoted as C[i]  (i being some integer) and are arbitray except that a suffcient number of  derivatives must be continuous.

deq = y D[u[x,y],x] + x D[u[x,y],y] == 0

so = DSolve[deq,u[x,y],{x,y}]

deq = y D[u[x,y],x] + x D[u[x,y],y] == 1

so = DSolve[deq,u[x,y],{x,y}]

deq = y D[u[x,y],x] - x D[u[x,y],y] == 0

so = DSolve[deq,u[x,y],{x,y}]

deq = x D[u[x,y,z],x] + y D[u[x,y,z],y] +  z D[u[x,y,z],z] == 0

so = DSolve[deq,u[x,y,z],{x,y,z}]

deq =  D[u[x,y],{x,2}] + D[u[x,y],{y,2}] == 0

so = DSolve[deq,u[x,y],{x,y}]

deq =  D[u[x,t],{x,2}] - 1/c^2 D[u[x,t],{t,2}] == 0

so = DSolve[deq,u[x,t],{x,t}]

PowerExpand[so]//Flatten

deq =  D[u[r,ϕ],{r,2}] + 1/r  D[u[r,ϕ],r] + 1/r^2 D[u[r,ϕ],{ϕ,2}]== 0

so = DSolve[deq,u[r,ϕ],{r,ϕ}]

11.3.2  Numerical Solution of Simple Partial Differential Equations

This solves the heat equation in one space dimension with the boundaryconditions  (Subscript[x, 0,1]) = 0, 5: T((Subscript[x, i]),t) = 0.
At the beginning a temperature peak is present around x= 2.5.

sol = NDSolve[{ D[T[x,t],t]==D[T[x,t],x,x]/10,T[x,0]==Exp[-(x-2.5)^2/0.05 ],T[0,t]==0,T[5,t]==0},T,{x,0,5},{t,0,4}]//Flatten

Plot3D[T[x,t]/.sol,{x,0,5},{t,0,4},AxesLabel-> {"x","t","T(x,t)  "},Ticks -> {Range[0,5],Range[0,4],Range[0,2]/2},
 PlotRange->All, BaseStyle -> FontSize -> 14]

11.3.3  Numerical Solution of Partial Differential Equations by 
	Difference or Finite Elements Methods

Numerical Solutions of Partial Differential Equations may be found by difference method or finite elements. 
For the first method, the partial derivatives are replaced by differences; in this way the differential equations are approximated by large linear systems of equations for the field values at the nodes of the mesh. In the book of 
Victor G. Ganzha and Evgenii V. Vorozhtsov:  Numerical Solutions for Partial Differential Equations:  
Problem Solving Using Mathematica, 
this is shown and programs are given. 

Prem K. Kythe, Pratap Puri, Michael R. Schaferkotter, Hiles K. Jones, Michael R. Schaeferkotter: 
Partial Differential Equations and Boundary Value Problems with Mathematica
CRC Press; 2 edition, 2002

11.4  	Exercises

11.1	A charged particle moves in static, homogeneous, crossed electric and magnetic fields and experiences a linear frictional 
	force.

d2r/dt2   =   q E  +  q (dr/dt  x  B)  - m a  dr/dt
      		
      		E = (0, E0, 0),  B = (0,0,B0),  a = const.

Initial data are: 	
		r(0)  = (0,0,0),    dr/dt(0)  =   (vx, vy, vz).

Solve the system of differential equations analytically (if possible) and numerically for suitable values of the parameters 
	and initial conditions.

11.2	Same problem as 1., but the frictional force is quadratic, i.e. it is

			- m a dr/dt | dr/dt |.
			
	Compute the trajectory by numeric solution for suitable values of the parameters and prepare a drawing.

11.3. 	A particle is projected from the ground with a given initial velocity  oblique to the field of gravity. Compute and draw the 
	trajectory. 
	Indicate equal time intervalls by dots on the trajectory. 
	Do the same for an additional linear or quadratic frictional force.

11.4.  	A 1-dimensional linear force   F  =  - x  acts on a  particle as long as   |x| < 1;  the force is zero for |x| > 1. 
	The particle starts with a given	initial speed at  x = 0. Compute and draw the position and the velocity as function of time 		and draw the phase space diagram for two initial velocities  such that the particle is bound and for two initial velocities 		belonging to  unbound motion.

11.5.  	Compute and draw the Poincare Map for the x,x'-plane of Sect. 11.2.6

11.6.   	Treat the following differential equation:	       y'' + 2 x y' + 2 y  =  0;	   
11.6.1 	Solve the differential equation analytically.
11.6.2 	Solve analytically the boundary values problem:	y(0) = 0, y(2) =  1. 
11.6.3  Solve numerically the initial value problem:	y(0) = 1, y'(0) = 0;   plot the solution.

11.7         Verify that   y[x] = Exp[x]     is  a solution of the following  differential equation:

   y[x]-x y'[x]+(-3+2 x) y''[x]+(2-x) Derivative[3][y][x] = 0.

11.8  	Compute the eigenvalue ν = 1 for the quantenmechanical Harmonic Oscillater by the 
	shooting method as it is done for ν = 0 in  subsubsection 11.2.7.1. Plot the curves for 
	⋁ around 1. Show the curvews of the shooting method in detail for  ν = 0.8, 0.9, 1.0 and 1.1 .