Maple: ODE and PDE system coupled

mscolli

Hello,

This is my first post and hopefully my question has not been answered elsewhere already as I realize it is annoying to answer the same type of posts over and over again.

I am working on a system of PDEs with one ODE, coupled. It is an SEIR model with one extra class for the group of people who have received treatment.

I believe I am having a problem because of the multiple dependencies in each of the equations and conditions.

Ideally I would like the system to be solved and plotted however I have no idea how to go about this with what I have here. I have searched the internet for other example but I have found nothing helpful. I have tried everything in the Maple help and used all options available to pdsolve. If somebody could take a look at the worksheet and guide me in the right direction that would be very helpful.

I will paste the contents of my worksheet below as well, I will attach it to this post if that is easier to read.

The version of Maple that I am using is Maple 10 student edition.

Thanks,
Shannon

> restart

> with(PDEtools); with(plots);

Warning, the name changecoords has been redefined

> with(DEtools)

> #Ideally I would like a system solved analytically and the results plotted. I would like the results for S,E,I,R,D plotted on the same graph.I have no idea what I have done wrong but I believe my issue lies in the fact that the system I am trying to solve is recursive.

>

>

> #alpha := 0.3257;

> #d := 0.00761;

> #mu := 0.0015;

> #omega := 1/30;

> #lambda := 0.01028;

> #gam := 0.49089;

> #R0 := 1.42;

> #p := 0.1;

> #beta := (R0 * d * (d + mu + gam) * (d + alpha))/(lambda*alpha);

> "#v:=0.8;"

> int(Ip(a, t), a = 0 .. infinity) := Ip(a, t);


> int(Del(a, t), a = 0 .. infinity) := Del(a, t);


> int(R(a, t), a = 0 .. infinity) := R(a, t);

int(Ip(a, t), a = 0 .. infinity) := Ip(a, t)

int(Del(a, t), a = 0 .. infinity) := Del(a, t)

int(R(a, t), a = 0 .. infinity) := R(a, t)

> "#ode:=diff(S(a,t),t) - lambda+d*S(a,t) + S(a,t)*beta*int(Ip(a,t),a=0..infinity)+ v*p*S(a,t) - omega*int(Del(a,t),a=0..infinity)-f*int(R(a,t),a=0..infinity)=0; "

> "#dsolve(ode, S(a,t));"

> [/ d \ / d \


> pdesys := [|--- E(a, t)| + |--- E(a, t)| + d E(a, t) + alpha E(a, t) = 0,


> [\ dt / \ da /


>


> / d \ / d \


> |--- Ip(a, t)| + |--- Ip(a, t)| + d Ip(a, t) + gam Ip(a, t) + mu Ip(a, t) = 0,


> \ dt / \ da /


>


> / d \ / d \


> |--- R(a, t)| + |--- R(a, t)| + f R(a, t) + d R(a, t) = 0,


> \ dt / \ da /


>


> / d \ / d \ ]


> |--- Del(a, t)| + |--- Del(a, t)| + d Del(a, t) = 0]


> \ dt / \ da / ]

> print();

[/ d \ / d \


[|--- E(a, t)| + |--- E(a, t)| + d E(a, t) + alpha E(a, t) = 0,


[\ dt / \ da /





/ d \ / d \


|--- Ip(a, t)| + |--- Ip(a, t)| + d Ip(a, t) + gam Ip(a, t) + mu Ip(a, t) = 0,


\ dt / \ da /





/ d \ / d \


|--- R(a, t)| + |--- R(a, t)| + f R(a, t) + d R(a, t) = 0,


\ dt / \ da /





/ d \ / d \ ]


|--- Del(a, t)| + |--- Del(a, t)| + d Del(a, t) = 0]


\ dt / \ da / ]

> IBCs := {E(a, 0) = 0, Ip(a, 0) = 0.05, R(a, 0) = 0, Del(a, 0) = 0,


>


> E(0, t) = beta Ip(0, t), Ip(0, t) = alpha E(0, t), R(0, t) = gam Ip(0, t),


>


> Del(0, t) = v p}

> print();

{E(0, t) = beta Ip(0, t), Ip(0, t) = alpha E(0, t), R(0, t) = gam Ip(0, t),





Del(0, t) = v p, E(a, 0) = 0, Ip(a, 0) = 0.05, R(a, 0) = 0, Del(a, 0) = 0}

> pdesol := pdsolve(pdesys)

> print();

{Del(a, t) = _F1(t - a) exp(-d a), E(a, t) = _F2(t - a) exp(-(d + alpha) a),





R(a, t) = _F3(t - a) exp(-(f + d) a),





Ip(a, t) = _F4(t - a) exp(-(d + gam + mu) a)}



>

Attached Files

preexposure.mw (105.9 KB, 2 views)