BVP4c für unbekannte Grenze lösen

Question

Ich versuche bvp4c zu verwenden, um ein System von 4 odes zu lösen. Das Problem ist, dass eine der Grenzen unbekannt ist.

Kann bvp4c damit umgehen? In meinem Code L ist das Unbekannte, für das ich mich löse.

Ich erhalte eine Fehlermeldung, die unten gedruckt wird.

function mat4bvp 
L = 8; 
solinit = bvpinit(linspace(0,L,100),@mat4init); 
sol = bvp4c(@mat4ode,@mat4bc,solinit); 
sint = linspace(0,L); 
Sxint = deval(sol,sint); 
end 

% ------------------------------------------------------------ 
function dtdpdxdy = mat4ode(s,y,L) 
Lambda = 0.3536; 
dtdpdxdy = [y(2) 
      -sin(y(1)) + Lambda*(L-s)*cos(y(1)) 
      cos(y(1)) 
      sin(y(1))]; 
end 
% ------------------------------------------------------------ 
function res = mat4bc(ya,yb,L) 
res = [ ya(1) 
     ya(2) 
     ya(3) 
     ya(4) 
     yb(1)]; 
end 
% ------------------------------------------------------------ 
function yinit = mat4init(s) 
yinit = [ cos(s) 
    0 
    0 
    0 
     ]; 
end 

Leider bekomme ich die folgende Fehlermeldung;

>> mat4bvp 
Not enough input arguments. 

Error in mat4bvp>mat4ode (line 13) 
      -sin(y(1)) + Lambda*(L-s)*cos(y(1)) 

Error in bvparguments (line 105) 
    testODE = ode(x1,y1,odeExtras{:}); 

Error in bvp4c (line 130) 
    bvparguments(solver_name,ode,bc,solinit,options,varargin); 

Error in mat4bvp (line 4) 
sol = bvp4c(@mat4ode,@mat4bc,solinit); 

Answer 1

Ein Trick einen variablen Endpunkt in eine festen ein zu transformieren ist, um die Zeitskala zu ändern. Wenn x'(t)=f(t,x(t)) die Differentialgleichung ist, gesetzt t=L*s, s aus 0 zu 1 und berechnet die zugehörigen Differentialgleichung für y(s)=x(L*s)

y'(s)=L*x'(L*s)=L*f(L*s,y(s)) 

nächsten Stich zu verwenden ist, um die globale Variable in einen Teil der Differentialgleichung zu transformieren, indem Berechnen Sie es als konstante Funktion. So ist das neue System

[ y'(s), L'(s) ] = [ L(s)*f(L(s)*s,y(s)), 0 ] 

und der Wert von L tritt als zusätzlichen Grenzwert.

---

Ich habe keine Matlab leicht verfügbar, in Python SciPy dies als

implementiert werden kann

from math import sin, cos 
import numpy as np 
from scipy.integrate import solve_bvp, odeint 
import matplotlib.pyplot as plt 

# The original function with the interval length as parameter 

def fun0(t, y, L): 
    Lambda = 0.3536; 
    #print t,y,L 
    return np.array([ y[1], -np.sin(y[0]) + Lambda*(L-t)*np.cos(y[0]), np.cos(y[0]), np.sin(y[0]) ]); 

# Wrapper function to apply both tricks to transform variable interval length to a fixed interval. 

def fun1(s,y): 
    L = y[-1]; 
    dydt = np.zeros_like(y); 
    dydt[:-1] = L*fun0(L*s, y[:-1], L); 
    return dydt; 

# Implement evaluation of the boundary condition residuals: 

def bc(ya, yb): 
    return [ ya[0],ya[1], ya[2], ya[3], yb[0] ]; 

# Define the initial mesh with 5 nodes: 

x = np.linspace(0, 1, 3) 

# This problem has multiple solutions. Try two initial guesses. 

L_a=8 
L_b=9 

y_a = odeint(lambda y,t: fun1(t,y), [0,0,0,0,L_a], x) 
y_b = odeint(lambda y,t: fun1(t,y), [0,0,0,0,L_b], x) 

# Now we are ready to run the solver. 

res_a = solve_bvp(fun1, bc, x, np.array(y_a.transpose())) 
res_b = solve_bvp(fun1, bc, x, np.array(y_b.transpose())) 

L_a = res_a.sol(0)[-1] 
L_b = res_b.sol(0)[-1] 
print "L_a=%.8f, L_b=%.8f" % (L_a,L_b) 

# Plot the two found solutions. The solution are in a spline form, use this to produce a smooth plot. 

x_plot = np.linspace(0, 1, 100) 
y_plot_a = res_a.sol(x_plot)[0] 
y_plot_b = res_b.sol(x_plot)[0] 

plt.plot(L_a*x_plot, y_plot_a, label='y_a, L=%.8f'%L_a) 
plt.plot(L_b*x_plot, y_plot_b, label='y_b, L=%.8f'%L_b) 
plt.legend() 
plt.xlabel("t") 
plt.ylabel("y") 
plt.show()