Re: Tricky differential equation
- To: mathgroup at smc.vnet.net
- Subject: [mg41553] Re: Tricky differential equation
- From: "Dr. Wolfgang Hintze" <weh at snafu.de>
- Date: Mon, 26 May 2003 05:46:30 -0400 (EDT)
- References: <bafqej$6rq$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Luiz, I have tried to find a appoximate solutions to your differential equation in terms of polynomials of degree m for the special boundary conditions u[0] = 0, u'[0] != 0. The method can be generalized. This is probably not a very elegant method in Mathematica but it seems to work. Perhaps someone in the group can improve the method. Here it is (notice that in the text below comments are not depicted explicitly; that's because I copied it from a notebook in which I used Ctrl+7 format for comments) In[362]:= ClearAll["Global`*"] Differential equation In[363]:= w[x_] := x*x u''[x] + x u'[x] + (x*x - 1) Sin[u[x]] Series ansatz In[364]:= u[x_] := a[0] + Sum[a[k] x^k, {k, 1, m}] Constructing a solution (e.g. for m=5; you can try other values of m and look what happens) by substituting the ansatz for u[x] into the differential eqaution and requiring that the coefficients of all powers of x vanish In[365]:= m = 5; In[366]:= cond = {}; sol = {}; For [i = 0, i <= m, AppendTo[cond, (D[w[t], {t, i}] /. t -> 0) == 0]; AppendTo[sol, a[i]]; i++]; You can display cond and sol if you like deleting the comment marks below (* cond sol *) In[370]:= asl = Solve[cond, sol] From In[370]:= Solve::"ifun": "Inverse functions are being used by \!\(Solve\), so some \ solutions may not be found." From In[370]:= Solve::"svars": "Equations may not give solutions for all \"solve\" \ variables." Out[370]= \!\({{a[5] -> \(a[1]\ \((20 + 40\ a[1]\^2 + 3\ a[1]\^4)\)\)\/3840, a[4] -> 0, a[3] -> \(-\(1\/48\)\)\ a[1]\ \((6 + a[1]\^2)\), a[2] -> 0, a[0] -> 0}, {a[5] -> 0, a[4] -> 0, a[3] -> 0, a[2] -> 0, a[1] -> 0, a[0] -> 0}}\) Substituting the coefficients into the ansatz In[371]:= h[x_] := u[x] /. asl Displaying the solution In[372]:= Clear[a] In[359]:= h[x] Out[359]= \!\({x\ a[1] - 1\/48\ x\^3\ a[ 1]\ \((6 + a[1]\^2)\) + \(x\^5\ a[1]\ \((20 + 40\ a[1]\^2 + 3\ \ a[1]\^4)\)\)\/3840, 0}\) Plotting the solution for a special value of a[1]: a[1] = 0.5; Plot[h[x][[1]], {x, 0, \[Pi]}, PlotRange -> {{0, \[Pi]}, {-1, 1}}]; or, for different a[1] For[z = 0.1, z < \[Pi], a[1] = z; Plot[h[x][[1]]/\[Pi], {x, 0, 2\[Pi]}, PlotRange -> {{0, \[Pi]}, {-1, 1}}]; z = z + 0.2]; Hope this hepls. Wolfgang Luiz Melo wrote: > Hello everyone, > > I'm trying to find the numerical solution of the following > differential equation (r is the independent variable): > > x''[r] + 1/r x'[r] + (p - 1/r^2)*Sin[x[r]]*Cos[x[r]] == 0 , > > with boundary conditions: x'[1] == 0 , and x[0] -> "has to be finite", > > but I'm having at least two problems: > > 1) I don't know how to submit the BC "finite" to Mathematica; > 2) The coefficient p is about 10^4. For this reason, it seems > that the Runge-Kutta method usually used for numerical > integration of ordinary differential equations turns out > to be unsuccessfull in our case. Do we need a special method > to solve this? > > The solution of this equation gives the internal magnetic structure > of a cylinder. The function x[r] is the angle between the > magnetization and the axial direction, and it depends on the radial > direction, r. > > I would like to plot the Cossine of the result as a function of r > (which varies from 0 to 1), for several values of p. > > Any help will be very appreciated! > Thank you > > Luiz Melo > > Ecole Polytechnique de Montreal, > Montreal, Quebec > luiz.melo at polymtl.ca > > > >