BVP4c解决未知边界问题

Question

我正在尝试使用bvp4c解决4个系统的问题。问题是其中一个边界是未知的。

bvp4c可以处理这个吗？在我的代码中，L是我正在解决的未知问题。

我在下面打印出错误信息。

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 

不幸的是我收到以下错误信息;

>> 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

将变量终点转换为固定终点的一个诀窍是更改时间刻度。如果x'(t)=f(t,x(t))是差分方程，设定t=L*s，s从0到1，并计算用于y(s)=x(L*s)

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

下一招采用相关的微分方程是全局变量转换成由差分方程的一部分将其计算为常量函数。因此，新的系统是

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

和L值出现额外的边界值。

---

我没有Matlab的一应俱全，在Python SciPy的这个可以作为

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()