你觉得有挑战性么【打靶法求初值并解微分方程】

356days

问题见附件图片。



我编的程序如下：

function  myfun()

cguess=0; % a guess at the value of y'(0)

actualvalue=fzero(@shooting,cguess)
% use ode45 to get correct solution using actualvalue and plot
tspan=[0 1];
y0=[0 actualvalue];
[h,y]=ode45(@concentrationprofile,tspan,y0);
plot(h,y(:,1),'r')
print -deps shootingexampleplot

function f=shooting(c)
% set up as a function, the input (c) is the initial slope y'(0)
% and the output (f) is the difference between the value of y at x=1
% and the number 1 (since y(1)=1 for the correct solution.
% This function is called by MATLABs fzero

% shooting.m
% trying to solve y''+ b1(h) * y' = 0 y(0)=0 y(1)=1
% using y(0)=0 and y’(0)=c
tspan=[0 1];
y0=[0 c];
[t,y]=ode45(@concentrationprofile,tspan,y0);
f=y(length(y),1)-1; % set function value to y(1)-1
return

function f=concentrationprofile(h,y)
% shootingexamplef.m
% the DE y''+ b1(h) * y' = 0 written as a system
q=quad(@(h) sin(h).^0.5,0,h);
s=sin(h)^0.5/q^(1/3);
b1=diff(s,h)./s

f(1)=y(2);
f(2)=-b1*y(2);
f=f(:);
return

当b1取常数时，程序没问题，但是当b1为附件图片所示时，就有问题了。

请教各位高手如何解决？