aylin

Thank you very much. Here, x[1], x[2], x[3] denote the concentration of uninfected hepatocytes, infected
hepatocytes, and free virus, respectively and they are should be positive. u1 and u2 are not parameters.The control functions u1(t) and u2(t) are bounded Lebesgue integrable functions.
The control u1(t) represents the eficiency of drug therapy in blocking
new infection, and the control u2(t) represents the eficiency of drug therapy
in inhibiting viral production.The control set defined by

{(u1,u2)|u1 and u2 are measurable, 0≤u1≤1 and 0≤u2≤1}

@ The problem is to maximize...

Answered: aylin 25

December 05 2014

0 0

@

The problem is to maximize the objective functional

J(u1(t); u2(t)) =int(x[1]-A[1]*(u[1])^2/2-A[2]*(u[2])^2/2,t=0..t_f)

where the parameters A[1]≥ 0 and A[2]≥ 0 are based on the benets and
costs of the treatment. Our target is to maximize the objective functional
by increasing the number of the uninfected cells, decreasing
the viral load, and minimizing the cost of treatment.

@ Thanks alot. The original contro...

Answered: aylin 25

December 05 2014

0 0

@

Thanks alot. The original control problem is in the following

dH/dt=λ-µ H-(1-u1)β H V+δ I,

dI/dt=(1-u1)β H V-σ I,

‎dV/dt=(1-u2)k I-γ‎V,

Please see the attached maple code.

I've used Pontryagin's maximum principle.

>	restart: unprotect('gamma');

>	L:=x[1]-A[1](u[1])^2/2-A[2](u[2])^2/2;

(1)

>	H:=L+psi[1](t)(lambda-mux[1]-(1-u[1])betax[1]x[3]+deltax[2])+psi[2](t)((1-u[1])betax[1]x[3]-sigmax[2]) +psi[3](t)((1-u[2])kx[2]-gamma*x[3]);

H := x[1]-(1/2)*(A[1]*u[1]^2)-(1/2)*(A[2]*u[2]^2)+psi[1](t)*(lambda-mu*x[1]-(1-u[1])*beta*x[1]*x[3]+delta*x[2])+psi[2](t)*((1-u[1])*beta*x[1]*x[3]-sigma*x[2])+psi[3](t)*((1-u[2])*k*x[2]-gamma*x[3])

(2)

>	du1:=diff(H,u[1]);

(3)

>	du2:=diff(H,u[2]);

(4)

>	ddu1:=-A[1]u[1]+psi[1](t)betax[1]x[3]-psi[2](t)betax[1]*x[3];

(5)

>	ddu2:=-A[2]u[2]-psi[3](t)k*x[2];

(6)

>	sol_u1 := solve(ddu1, u[1]);

sol_u1 := beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1]

(7)

>	sol_u2 := solve(ddu2, u[2]);

sol_u2 := -psi[3](t)*k*x[2]/A[2]

(8)

>	Dx2:=subs(u[1]=betax[1]x[3](psi[1](t)-psi[2](t))/A[1],u[2]=-psi[3](t)k*x[2]/A[2], H );

Dx2 := x[1]-(1/2)*(beta^2*x[1]^2*x[3]^2*(psi[1](t)-psi[2](t))^2/A[1])-(1/2)*(psi[3](t)^2*k^2*x[2]^2/A[2])+psi[1](t)*(lambda-mu*x[1]-(1-beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1])*beta*x[1]*x[3]+delta*x[2])+psi[2](t)*((1-beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1])*beta*x[1]*x[3]-sigma*x[2])+psi[3](t)*((1+psi[3](t)*k*x[2]/A[2])*k*x[2]-gamma*x[3])

(9)

>	ode1:=diff(psi[1](t),t)=-diff(H,x[1]);

(10)

>	ode2:=diff(psi[2](t),t)=-diff(H,x[2]);

(11)

>	ode3:=diff(psi[3](t),t)=-diff(H,x[3]);

(12)

>	restart: #Digits:=10: unprotect('gamma'); lambda:=510^5: mu:=0.003: beta:=410^(-10): delta:=0: alpha:=0.043: sigma:=alpha+delta: k:=6.24: gamma:=0.65: A[1]:=10: A[2]:=2: u[1]:=betax[1](t)x[3](t)(psi[1](t)-psi[2](t))/A[1]: u[2]:=-psi[3](t)k*x[2](t)/A[2]:

>	ics := x[1](0)=1.7*10^8, x[2](0)=0,x[3](0)=400,psi[1](20)=0,psi[2](20)=0,psi[3](20)=0:

>

ode1:=diff(x[1](t), t)=lambda-mu*x[1](t)-(1-u[1])*beta*x[1](t)*x[3](t)+delta*x[2](t),
diff(x[2](t), t) =(1-u[1])*beta*x[1](t)*x[3](t)-sigma*x[2](t),
diff(x[3](t), t) =(1-u[2])*k*x[2](t)-gamma*x[3](t),
diff(psi[1](t), t) =-1+mu*psi[1](t)+beta*x[3](t)*(1-u[1])*(psi[1](t)-psi[2](t)),
diff(psi[2](t), t) =delta*psi[1](t)+sigma*psi[2](t)-psi[3](t)*(1-u[2])*k,
diff(psi[3](t), t) = beta*x[1](t)*(psi[1](t)-psi[2](t))*(1-u[1])+gamma*psi[3](t);

(13)

>	sol:=dsolve([ode1,ics],numeric, method = bvp[midrich]);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

>

Download mapleprime.mws

transversality conditions...

Answered: aylin 25

October 23 2014

0 1

@Carl Love

λ1,λ2,λ3 are transversality conditions where

λi(t_end)=0

and t_end=70.

full code...

Answered: aylin 25

July 13 2014

0 0

Thanks for your response

my full code is following.

restart:

unprotect('gamma');
x[0]:=140: y[0]:=145:v[0]:=18:
A1:= 900: A2:=1000:
mu:=0.73:
gamma:=250:
beta:=0.0014:
r:=0.01:
delta:=0:
a:=0.0693:
k:=340:

#Step 1
n:= 100:
lambda__1[n]:= 0: lambda__2[n]:= 0: lambda__3[n]:= 0:
u__1[0]:= 0: u__2[0]:= 0:

#Step 2:
for i from 0 to n-1 do
x[i+1]:= x[i]+h*(r*x[i+1]*(1-(x[i+1]+y[i])/k)-(1-u__1[i])*beta*v[i]*(x[i+1]/(x[i+1]+y[i])));
y[i+1]:= y[i]+h*(((1-u__1[i])*beta*v[i]*x[i+1]/(x[i+1]+y[i+1]))-a*y[i+1]);
v[i+1]:= v[i]+h*((1-u__2[i])*gamma*y[i+1]-mu*v[i+1]);
lambda__1[n-i-1]:=lambda__1[n-i]+h*(-1-lambda__1[n-i-1]*(r*(1-(x[i+1]+y[i])/k)
-r*x[i+1]/k-(1-u__1[i])*beta*v[i+1]*y[i+1]/(x[i+1]+y[i])^2)-lambda__2[n-i]*(1-u__1[i])*beta*v[i+1]*y[i+1]/(x[i+1]+y[i])^2);
lambda__2[n-i-1]:= lambda__2[n-i]+h*(lambda__1[n-i-1]*(r*x[i+1]/k-(1-u__1[i])*beta*v[i+1]*x[i+1]/(x[i+1]+y[i])^2)-lambda__2[n-i-1]*((1-u__1[i])*beta*v[i+1]*x[i+1]/(x[i+1]+y[i])^2+a)-lambda__3[n-i]*gamma*(1-u__2[i]));
lambda__3[n-i-1]:= lambda__3[n-i]+h*((1-u__1[i])*beta*x[i+1]/(x[i+1]+y[i])*(lambda__1[n-i-1]-lambda__2[n-i-1])+mu*lambda__3[n-i-1]);
R1:= (1/A1*(x[i+1]+y[i+1]))*(lambda__1[n-i-1]-lambda__2[n-i-1])*beta*v[i+1]*x[i+1];
R2:= -(1/A2)*lambda__3[n-i-1]*gamma*y[i+1];
u__1[i+1]:= min(1, max(R1,0));
u__2[i+1]:= min(1, max(R2,0));
end do:

When I run this code in maple I am facing with "Error, recursive assignment".

Please help me to solve this equations on Maple.