Code covered by the BSD License

# Chebpack

### Damian Trif (view profile)

15 Jul 2011 (Updated )

The MATLAB package Chebpack solves specific problems for differential or integral equations.

[t1opt,I1opt,F]=ex3
```function [t1opt,I1opt,F]=ex3
% Gompertz equation x'=x(r-gamma*log(x)),t in [0,5]
%                   x(0)=0.2, r=0.03, gamma=1
% Find t1 and I1 such that after the removal biomass x(t1+0)=x(t1)-I1,
% F(t1,I1)=I1+x(5) be maximal
% From A. Dishliev et. al., Specific asymptotic properties of the solutions
% of impulsive differential equations. Methods and applications, Academic
% Publications, Ltd., 2012
% use: [t1,I1,F]=ex3;
%
r=0.03;gamma=1;x0=0.2;tfin=5;n=32;F=zeros(20,10);
for i=1:20
t1(i)=2+i*0.02;
[~,x]=gompertz(0,t1(i),x0);a1(i)=real(x(end));
for j=1:10
I1(j)=0.75+j*0.01;if I1(j)>=a1(i), break;end
[~,x]=gompertz(t1(i),tfin,a1(i)-I1(j));a2(i,j)=real(x(end));
F(i,j)=I1(j)+a2(i,j);
end
end
[X,Y]=meshgrid(t1,I1);
figure(1);[C,H]=contour(X,Y,F',[1.655,1.66,1.664,1.666]);hold on;
xlabel('t_1');ylabel('I_1');axis([2 2.4 0.75 0.85]);grid;
clabel(C,H);
plot(2.2012,0.8072,'*r');
[X,F]=fminsearch(@gpzopt,[2,0.75]);
t1opt=X(1);I1opt=X(2);F=-F;
figure(2);
[tl,xl]=gompertz(0,t1opt,x0);plot(tl,xl);hold on;a1=real(xl(end));
[tr,xr]=gompertz(t1opt,tfin,a1-I1opt);plot(tr,xr);hold on;
xlabel('t');ylabel('x(t)');grid;
% Calculate the residual for the optimal solution
Xl=[x2t(xl,kind);zeros(n,1)];Xr=[x2t(xr,kind);zeros(n,1)];
vl=t2x(Xl,kind);vr=t2x(Xr,kind);
Dl=deriv(2*n,[0,t1opt]);Dr=deriv(2*n,[t1opt,tfin]);
resl=t2x(Dl*Xl-r*Xl,kind)+gamma*vl.*log(vl);
resr=t2x(Dr*Xr-r*Xr,kind)+gamma*vr.*log(vr);
figure(3);
semilogy([pd(2*n,[0,t1opt],kind);pd(2*n,[t1opt,tfin],kind)],abs([resl;resr]));
grid;xlabel('t');ylabel('residual')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F=gpzopt(X)
tt=X(1);II=X(2);
[~,xl]=gompertz(0,tt,x0);aa1=real(xl(end));
[~,xr]=gompertz(tt,tfin,aa1-II);aa2=real(xr(end));
F=-(II+aa2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t,x]=gompertz(t0,t1,x0)
dom=[t0,t1];kind=2;t=pd(n,dom,kind);D=deriv(n,dom);T=cpv(n,dom(1),dom);
A=D-r*speye(n);A(n,:)=T;
w=zeros(n,1);w(1)=-gamma*x0*log(x0);
b=x2t(w,kind);b(end)=x0;
xold=A\b;
for k=1:50
xval=t2x(xold,kind);w=-gamma*xval.*log(xval);
b=x2t(w,kind);b(end)=x0;xnew=A\b;
var=norm(xnew-xold);if var<1.e-12,break;end
xold=xnew;
end
x=xval;
end
end```