1
$\begingroup$

I am trying to solve below 3rd order DE using Runge Kutta methed but we have problem at t=0,

ClassicalRungeKuttaCoefficients[4, prec_] := 
  With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
     bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
   N[{amat, bvec, cvec}, prec]]

σ = .07;
ρ = 1000; g = 10; σT = .3*10^-3; 
lσ = Sqrt[σ/(ρ*g)]; ΔT = .02;
Cr = (σT*ΔT)/σ; h0 = 1*10^-3;
Amp = (3*lσ^2*Cr)/(2*h0^2); T0 = 300;
T = ΔT*Exp[-(t/(200*10^-6))^2];
θ = (T - T0)/ΔT;
solution1 = 
   NDSolve[{y''[t] + 1/t*y'[t] - y[t] == 
     Amp*Integrate[θ/(y[t] + 1), {t, -.008, .008}], 
      y[-0.008] == 0, y[0.008] == 0}, y[t], {t, -.008, .008}, 
     Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
      "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
     StartingStepSize -> 1/10000]

these error appeared, and plot is not symmetric around t=0

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ Not sure I understand the question. Your equation contains the term y'[t]/t, so I'm not really sure why it's surprising to find an error at t=0. But the solution1 you get out seems pretty useable as far as I can see. What exactly is the problem at t=0 (since you get a solution at the end)? What are you wanting to do with the solution that you can't do now? $\endgroup$
    – aardvark2012
    Jul 22, 2017 at 11:44
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Jul 23, 2017 at 4:15
  • $\begingroup$ It looks like there's a contradiction between the boundary conditions at the top of your post, $h(\pm\infty)=0$, and the solution graph where $h(r)$ tends towards 1. Should the BCs be $h'(\pm\infty)=0$? $\endgroup$
    – Chris K
    Jul 29, 2017 at 13:34

1 Answer 1

Reset to default
1
$\begingroup$

Because the ODE is singular at t == 0, it seems prudent to solve the ODE with DSolve rather than NDSolve. Unfortunately, using DSolve directly fails. So, instead solve

s = DSolveValue[{y''[t] + 1/t*y'[t] - y[t] == c, 
      y[-.008] == 0, y[.008] == 0}, y[t], {t, -.008, .008}] // Chop
(* c (-1 + 0.999984 BesselJ[0, I t]) *)

where c is the value of the integral, determined as follows.

ss = s/c;
int[c_] := Amp NIntegrate[θ/(ss c + 1), {t, -.008, .008}]
sc = c /. FindRoot[int[c] == c, {c, 1}, Evaluated -> False] // Chop
(* -0.211214 *)

Thus the solution is ss sc

Plot[sc ss, {t, -.008, .008}, ImageSize -> Large, 
    AxesLabel -> {t, y}, LabelStyle -> Directive[Bold, Black, 12]]

enter image description here

Improve this answer
$\endgroup$
11
  • $\begingroup$ @MichaelE2 I share your discomfort with the answer. However, the only way to eliminate the discontinuity is to force the coefficient of BesselY to be zero, in which case the given boundary conditions cannot be satisfied. $\endgroup$
    – bbgodfrey
    Jul 22, 2017 at 19:53
  • $\begingroup$ Thanks , typo in range corrected(-..008 to .008), But after that y(t) does not goes to Asymptotically zero.We need only real solution. We can change boundary condition 0 to some finite no. $\endgroup$
    – Gopal Verma
    Jul 22, 2017 at 20:35
  • $\begingroup$ @GopalVerma That change makes a big difference. $\endgroup$
    – bbgodfrey
    Jul 22, 2017 at 20:48
  • $\begingroup$ But,It should goes to Asymptotically zero. $\endgroup$
    – Gopal Verma
    Jul 22, 2017 at 21:32
  • 1
    $\begingroup$ I guess my previous comment is academic now. I like the sc = FixedPoint[int, -1.] trick. BTW, you can get rid of the last imaginary I with FunctionExpand[s]. $\endgroup$
    – Michael E2
    Jul 22, 2017 at 23:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.