function prob_3_26 % % Written by: % -- % John L. Weatherwax 2005-05-07 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- global c alpha beta infty options = []; options = bvpset(options,'Vectorized','on'); options = bvpset(options,'FJacobian',@odeJac); options = bvpset(options,'BCJacobian',@bcJac); color = [ 'r', 'b', 'g', 'm', 'k' ]; wave_speed = [5, 10, 15, 20, 25]; for i = 1:5 c = wave_speed(i); alpha = (-c + sqrt(c^2 + 4))/2; beta = (-c + sqrt(c^2 - 4))/2; infty = 10*c; % if i == 1 solinit = bvpinit(linspace(-infty,infty,20),@guess); % else % solinit = bvpinit(sol,[-infty,infty]); % end sol = bvp4c(@ode,@bc,solinit,options); figure; plot(sol.x,sol.y(1,:),color(i)); hold on ; gus = guess(sol.x); plot(sol.x,gus(1,:),'xk'); axis([-250 250 0 1]); drawnow end %legend('c = 5', 'c = 10', 'c = 15', 'c = 20', 'c = 25',3); %print -depsc ch3fig5 %====================================================== function v = guess(z) global c alpha beta infty METHOD=2; switch METHOD case 1 if z > 0 v = [exp(beta*z); beta*exp(beta*z)]; else v = [(1 - exp(alpha*z)); -alpha*exp(alpha*z)]; end case 2 tmp = 1+exp(z/c); v = [ 1 ./ tmp; -(1/c) * (tmp-1) ./ (tmp.^2); ]; otherwise end function dydz = ode(z,y) global c alpha beta infty dydz = [ y(2,:); -(c*y(2,:) + y(1,:).*(1 - y(1,:))) ]; function dFdy = odeJac(z,y) global c alpha beta infty dFdy = [ 0, 1 (-1 + 2*y(1)), - c ]; function res = bc(ya,yb) global c alpha beta infty res = [ ya(2)/(ya(1) - 1) - alpha yb(1)/exp(beta*infty) - 1 ]; % deriveFbc.m derives these partial derivatives. function [dBCdya, dBCdyb] = bcJac(ya,yb) global c alpha beta infty dBCdya = [ -ya(2)/(ya(1) - 1)^2, 1/(ya(1) - 1) 0 0 ]; dBCdyb = [ 0 0 1/exp(beta*infty), 0 ];