lsq.txt

打靶法求微分方程边值问题

%打靶法求微分方程边值问题
% y1'=y2(x),y2'=(p*hc/A*k)*(y1(x)-293)
% lsq
while 1
y2=input('Type gradient,y2(0);or -99999 to quit:');
if y2<-88888 break, end
A=0.0001;p=0.01;hc=120;k=60;b=293;
a=p*hc/A/k;
n=1;x(1)=0;h=0.01;
y(:,1)=[493;y2];
while x<=0.3
k1=h*f_shoot(y(:,n),x(n),a,b);
k2=h*f_shoot(y(:,n)+k1/2,x(n)+h/2,a,b);
k3=h*f_shoot(y(:,n)+k2/2,x(n)+h/2,a,b);
k4=h*f_shoot(y(:,n)+k3,x(n)+h,a,b);
y(:,n+1)=y(:,n)+(1/6)*(k1+2*k2+2*k3+k4);
x(n+1)=n*h;
if(x(n)-0.2001)*(x(n)-0.1999)<0
y2_end=y(2,n+1),break
end
n=n+1;
end
end
plot(x,y(1,,'--',x,y(2,/10,':')
function f=f_shoot(y,x,a,b)
f=[y(2);a*(y(1)-b)]; 

end
plot(t,y(1:3,);