0

I have solved the following bvp problem using bvp solver in python.

d4y/dx4= 0.00033*V/(0.000001-y)^(2) , y(0)=y'(0)=y(1)=y'(1)=0 In the above eqn 'V' is a parameter which has been varied using the for loop. The interesting part is that the solution to the above differential equation should be unstable for V>Vo. The bvp solver still yields some wrong values for V>Vo. How do I make the solver stop computing as soon as this instability arises?

Anweshan
  • 29
  • 6
  • 1
    How do you measure "instability", theoretically and algorithmically? Do you evaluate the `.success` field of the solution structure? – Lutz Lehmann May 04 '20 at 08:33
  • No, I did not evaluate it yet using the code. I wanted to know, whether the stability of this differential equation could be checked in matlab or python. – Anweshan May 04 '20 at 08:40
  • What is "stability"? The BVP solvers transform your problem into a huge, sparse and slightly non-linear system of equations that is then solved with some non-linear solver. Then an error estimation for the local error tells if the grid needs to be refined, and if yes, the process is started again. This can fail unrecoverably at many points, leading to error messages like "singular Jacobian" (of this large system) or "maximum node number exceeded", which in most cases means that the current trial solution is far away from any actual solution. Is any of that related to your "stability"? – Lutz Lehmann May 04 '20 at 08:52
  • The equation has only stable solution for VVo, I should not receive any solution at all. I would even be happy to see that the solver could not converge to a proper solution. But what I am getting is some solution which is not expected at all – Anweshan May 04 '20 at 09:21
  • Again, did you check the `.success` field of the solver result structure? Print out the `.message` field? – Lutz Lehmann May 04 '20 at 09:29
  • Could you please a suggest a way to use the .message and .success field? I have the code for the above differential eqn compiled but don't have an idea about how the functions .success and .message can be called – Anweshan May 04 '20 at 09:50

1 Answers1

1

For the normalized equation (changed scale of y and V)

    y''''*(1e-6-y)**2 = 3.3e-4*V
    (1e6*y)''''*(1-1e6*y)**2 = 3.3e14*V

    u = 1e6*y,   c = 3.3e14*V

    u'''' = c/(1-u)**2

I get a critical value for c=70.099305, that is, V0=0.2124e-12. For very small c the solution is likewise small and close to y(t)=c/24*(t*(1-t))**2. For c close to the critical value the grid refinement concentrates at the maximum close to y=0.4.

c=70.099305

def f(t,u): return [u[1],u[2],u[3],c/(1-u[0])**2]
def bc(u0,u1): return [u0[0], u0[1], u1[0], u1[1]]

t = np.linspace(0,1,5);
u = np.zeros([4,len(t)])
res = solve_bvp(f,bc,t,u, tol=1e-4, max_nodes=5000)
print(res.message)

%matplotlib inline
if res.success:
    plt.figure(figsize=(10,5))
    t = np.linspace(0,1,502)
    plt.plot(t, c/24*(t*(1-t))**2,c='y', lw=3)
    plt.plot(t,res.sol(t)[0],'b')
    plt.plot(res.x,res.y[0],'xr')
    plt.grid(); plt.show()

plot of solution and small parameter approximation

blue - numerical solution, yellow - approximation for small c

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • The rescaled maximum of `y` is at `0.4e-6`, if the unit is [m], then this is `0.4 [um]`. That the number of required mesh nodes explodes by changing from a few hundreds for the the `c` used above to more than 5000 (or 50000) for a very slightly larger `c` indicates that there is no solution that is close to the given one. I still do not know what kind of "stability" you are exploring, so I can not say how the difficulties of the solver relate directly to your "stability". – Lutz Lehmann May 04 '20 at 11:55
  • 1
    Thank you for your valuable guidance. The code that you suggested is giving perfect results as expected. Actually I was making an error in boundary conditions. It was my bad completely. Once again, thank you. Your code works perfectly fine. – Anweshan May 04 '20 at 12:53