打靶法求边值的matlab实现

参考代码：

% 解的初始估计

solinit = bvpinit(linspace(0,1,100),[1 -50])

 

% BVP问题求解

ode = @(t,x) [x(2) -pi^2/4*(1-x(1))]

bc = @(ya,yb)[ya(2)-pi yb(2)+pi/2]

sol = 前简bvp4c(ode,bc,solinit)

 

% 结果绘图

t=sol.x

vars={'x', '搭悔誉x'''}

for i=1:length(vars)

    subplot(2,1,i)

    plot(t,sol.y(i,:))

    知段xlabel('t')

    ylabel(vars{i})

end

4个微分方程应该只有4个约束（包括初值或终值）吧？

但按照你给的初值和终值，一共有6个约束。

可以庆颤这样理解：如果4个初誉扒败值都给定了，那么微分方程就有充足的条件可以直接求解了，不一定保证终此如值符合你的要求。

第一个例子

导d仿真Matlab源代码(missile simulation)

http://www.matlabsky.com/thread-284-1-1.html

第二例子

%by dynamic

%see also http://www.matlabsky.com

%2009.2.18

%比例导引理想d道胡改仿真凳做备Matlab代码

%给出了比例导引法的差分方程, 采用Matlab语言, 对比例导引法的导d理想d枣毁道进行了三维数据仿真, 绘制出三维理想d道。随导d与目标参数变化情况, 对理想d道的特性进行了分析。

clear all

clc

clear

tt=0.1

sm=0.6*tt

st=0.42*tt

x(1)=0y(1)=0z(1)=0

pmr(:,1)=[x(1)y(1)z(1)]

ptr(:,1)=[25510]

m=3

q(1)=0

o(1)=0

a(1)=0

for(k=2:600)

ptr(:,k)=[25-0.42*cos(pi/6)*tt*k510+0.42*sin(pi/6)*k*tt]

r(k-1)=sqrt((ptr(1,k-1)-pmr(1,k-1))^2+(ptr(2,k-1)-pmr(2,k-1))^2+(ptr(3,k-1)-pmr(3,k-1))^2)

c=sqrt((ptr(1,k)-pmr(1,k-1))^2+(ptr(2,k)-pmr(2,k-1))^2+(ptr(3,k)-pmr(3,k-1))^2)

b=acos((r(k-1)^2+st^2-c^2)/(2*r(k-1)*st))

dq=acos((r(k-1)^2-st^2+c^2)/(2*r(k-1)*c))

if abs(imag(b))>0

b=0.0000001

end

if abs(imag(dq))>0

dq=0.0000001

end

q(k)=q(k-1)+dq

o(k)=o(k-1)+m*dq

a(k)=o(k)-q(k)

c1=r(k-1)*sin(b)/sin(a(k)+b)

c2=r(k-1)*sin(a(k))/sin(a(k)+b)

c3=sqrt((c1-sm)^2+(c2-st)^2+2*(c1-sm)*(c2-st)*cos(a(k)+b))

dq=a(k)-acos(((c1-sm)^2+c3^2-(c2-st)^2)/(2*(c1-sm)*c3))

if abs(imag(dq))>0

dq=0.0000001

end

q(k)=q(k-1)+dq

o(k)=o(k-1)+m*dq

a(k)=o(k)-q(k)

c1=r(k-1)*sin(b)/sin(a(k)+b)

c2=r(k-1)*sin(a(k))/sin(a(k)+b)

c3=sqrt((c1-sm)^2+(c2-st)^2+2*(c1-sm)*(c2-st)*cos(a(k)+b))

dq=a(k)-acos(((c1-sm)^2+c3^2-(c2-st)^2)/(2*(c1-sm)*c3))

if abs(imag(dq))>0

dq=0.0000001

end

q(k)=q(k-1)+dq

o(k)=o(k-1)+m*dq

a(k)=o(k)-q(k)

c1=r(k-1)*sin(b)/sin(a(k)+b)

c2=r(k-1)*sin(a(k))/sin(a(k)+b)

c3=sqrt((c1-sm)^2+(c2-st)^2+2*(c1-sm)*(c2-st)*cos(a(k)+b))

x1(k)=ptr(1,k-1)+c2/st*(ptr(1,k)-ptr(1,k-1))

y1(k)=ptr(2,k-1)+c2/st*(ptr(2,k)-ptr(2,k-1))

z1(k)=ptr(3,k-1)+c2/st*(ptr(3,k)-ptr(3,k-1))

x(k)=pmr(1,k-1)+sm/c1*(x1(k)-pmr(1,k-1))

y(k)=pmr(2,k-1)+sm/c1*(y1(k)-pmr(2,k-1))

z(k)=pmr(3,k-1)+sm/c1*(z1(k)-pmr(3,k-1))

pmr(:,k)=[x(k)y(k)z(k)]

r(k)=sqrt((ptr(1,k)-pmr(1,k))^2+(ptr(2,k)-pmr(2,k))^2+(ptr(3,k)-pmr(3,k))^2)

if r(k)<0.06

break

end

end

sprintf('遭遇时间：%3.1f',0.1*k)

figure(1)

plot3(pmr(1,1:k),pmr(2,1:k),pmr(3,1:k),'k',ptr(1,:),ptr(2,:),ptr(3,:))

axis([0 25 0 5 0 25])

text(x(80),y(80),z(80),'\leftarrow 比例导引')

grid on

---