1

I cannot find out where the difference comes from.

First pde, solBM30, has variable boundary condition which I evaluate at a certain point. The second pde, solBM32 is identical, except it has a fixed boundary condition which is then evaluated at the same point. The results should be equal right?

σ = 0.; 
K = 110; 
S = 20; 
E = 10; 
θ = 0.000023637253287033743; 
κ = 0.008622227075117428; 
cr = -0.00024031738413887727; 
small = 0.01; 

solBM30 = 
  NDSolve[
    {D[F[r0, t], {t}] == (-(θ - κ r0)) D[F[r0, t], {r0}], F[r0, K] == 10000 r0}, 
    F, {t, 0, K}, {r0, cr - small, cr + small}]; 

solBM32 = 
  NDSolve[
    {D[W[r0, t], {t}] == (-(θ - κ r0)) D[W[r0, t], {r0}], W[r0, K] == 10000 cr}, 
    W, {t, 0, K}, {r0, cr - small, cr + small}];

(*Plot[slice[r0],{r0,rd,ru}]*)
F[cr, 0] /. solBM30
W[cr, 0] /. solBM32

Edit: Added the 'full' PDE. Now I get another error.

solBM1 = NDSolve[{0 == 
    0.5*σ^2*D[V[r0, t], {r0}, {r0}] + (θ - κr0)*
      D[V[r0, t], {r0}] - r0*V[r0, t] + D[V[r0, t], {t}], 
   V[r0, K] == 
    1, (D[V[r0, t], {t}] + θ*D[V[r0, t], {r0}] /. r0 -> 0) == 
    0}, {V}, 
     {t, 0, K}, {r0, cr - small, cr + small}]

Edit2:

r1 = -0.000240317
M = 111;
si = 0.000212693;
b = 0.00862223;
a = 0.0000236373;
R = 1/10;

eq = 0 == 0.5*si^2*D[u[t, x], {x, 2}] + (a - b*x)*D[u[t, x], x] - 
    x* u[t, x] - D[u[t, x], t];
ic = u[0, x] == 1 ;
bc = {(-D[u[t, x], {t}] + a*D[u[t, x], {x}] /. x -> 0) == 0, 
  u[t, R] == Exp[-R*t]}
domain = {0, R};
difforder = 4;
points = 25;
grid = Array[# &, points, domain];
(*Definition of pdetoode isn't included in this post,please find it \
in the link above.*)
ptoofunc = pdetoode[u[t, x], t, grid, difforder];
removeredundant = #[[2 ;; -2]] &;
ode = removeredundant@ptoofunc@eq;
odebc = ptoofunc@bc
odeic = ptoofunc@ic;
tend = M;
sollst = NDSolveValue[{ode, odeic, odebc}, u /@ grid, {t, 0, tend}];
sol = rebuild[sollst, grid];

Plot3D[sol[t, x], {t, 0, tend}, {x, ##}] & @@ domain

Edit3 (Beware I corrected a mistake in the first b.c in edit2 ):

vas[t_, r0_] = a/b + (r0 - a/b)*Exp[-b*t];
B[t_, T_] := (1 - E^(-b (T - t)))/b
A[T_, t_] := (a/b - si^2/(2 b^2)) (B[t, T] - (T - t)) - (
  si^2 B[t, T]^2)/(4 b)
Z1[t_, T_] := E^(A[T, t] - 0*B[t, T])
Z2[t_, T_] := E^(A[T, t] - vas[t, 0] B[t, T])
Plot[{sol[M - t, 0], Z1[t, M], Z2[t, M]}, {t, 0, M}]
  • Do not use E, K - Built-in Symbols. – Alex Trounev Mar 07 '19 at 16:51
  • 1
    Here is the hyperbolic equation, the slope of the characteristic changes sign on the left and right border. Therefore, it is necessary to put two boundary conditions. – Alex Trounev Mar 07 '19 at 16:57
  • It's necessary to provide b.c. in r0 direction when solving the problem numerically. Related: https://mathematica.stackexchange.com/q/73961/1871 https://scicomp.stackexchange.com/q/28744/5331 – xzczd Mar 08 '19 at 06:54
  • Many thanks! This means that a boundary condition must be defined as a function of t? Can I choose any r0 to define this function ? Are there any other ways to find a b.c? – Øyvind Foshaug Mar 08 '19 at 09:30
  • I believe this article gives me guidance: link – Øyvind Foshaug Mar 08 '19 at 12:46
  • I edited the post to reflect updates based on the help you gave and the article. This is the pde I am aiming at but with variable boundary condition such as the stripped down version I first gave. But first Id like the b.c to work. – Øyvind Foshaug Mar 08 '19 at 13:09
  • You need to add @xzczd in your comment, or I won't get the reminder. Then, 1. The full PDE is completely different from the simplified one, they belong to 2 different categories. 2. One more b.c. approximating b.c. at infinity is necessary. 3. NDSolve cannot handle the b.c. at r0==0, so we need to discretize the equation all by ourselves, pdetoode can be used for the task. – xzczd Mar 08 '19 at 15:32
  • @xzczd. Does that mean that if 2. is solved, no discretization is needed (seems like a lengthy task)? – Øyvind Foshaug Mar 08 '19 at 18:03
  • @xzczd , pdetoode seems promising. Do you have an exmaple similar to the above equation? – Øyvind Foshaug Mar 08 '19 at 22:20
  • No, issue 2 and issue 3 are separate i.e. we need to add the missing boundary and discretize the equation ourselves. We already have several questions similar to yours in this site, for example this and this, you may have a look. – xzczd Mar 09 '19 at 06:09
  • @xzczd, please see edit2. I am having trouble using pdetoode. I am new to MMA so my problems may be trivial. Could you also explain the difforder of x (in the first example it is 4, so I recon it is 4 in my application - but I am unsure here). Otherwise the first example you gave seems to match my application of pdetoode. – Øyvind Foshaug Mar 23 '19 at 16:34
  • Use R=1/10 to define R, MachinePrecision number can cause trouble because pdetoode currently uses === to check if an equation is on the boundary. (Well, perhaps I should re-design this part. ) 2. domain should be domain={0,R}. 3. difforder is short for difference order i.e. it's the order of the difference formula used for discretization.
  • – xzczd Mar 24 '19 at 03:10
  • I've updated pdetoode, now it can handle definition like R=0.1 properly. – xzczd Mar 24 '19 at 03:29
  • @xzczd. Thanks alot! I am sorry if this seems like a moving target. My problem is now to match the analytical solution. In Edit3 A and B are the solutions to the Riccatti equations. When r0>0 r0 can be a constant (Z1) and the analytical and numerical yields the same. But when r <= 0 the correct solution is only achieved with r as a function (Z2). Now the funny part. The numerical and the anaytical (Z1) are identical and seems to commit the same mistake. Z2 is what you expect - can be shown with a simple ODE. It would be interesting to see if you see what is going on here. Thank you! – Øyvind Foshaug Mar 24 '19 at 18:18
  • I'm sorry, but I don't understand the question. Do you mean you've deduced 2 analytic solutions Z1 and Z2, but Z1 is wrong and the numeric solution is the same as Z1? If so, the only thing I know at this point is, I believe the numeric solution is reliable enough, if it doesn't match the analytic solution, you should check whether same problem is solved or not. – xzczd Mar 25 '19 at 05:27
  • @xzczd. I see. Solution Z2 is correct, The numerical solution is not. By accident I achieve the wrong numerical solution by using a constant r=r(0) in B (some references suggested this). I was hoping this could be a clue as to why the numerical solution is off. The problem I believe has nothing to do with Pdetoode. Normal NDSolve also fives this result. Pdetoode provided nice help in supplying more b.c's, but to no avail. Let me know if I can provide you with more information. – Øyvind Foshaug Mar 25 '19 at 10:09
  • @xzczd. There is one more b.c. to try: D[u[t, x], {t}] /. x -> 0==0. Pdetoode then yields an error. Could you please have a look?Thank you ! – Øyvind Foshaug Mar 25 '19 at 10:13
  • It should be (D[u[t, x], {t}] /. x -> 0)==0 or D[u[t, x], {t}] ==0/. x -> 0. – xzczd Mar 25 '19 at 16:42
  • Thanks @xzczd. Still I dont get it work I am afraid. – Øyvind Foshaug Mar 27 '19 at 07:05