Maple: ODE and PDE system coupled

  1. 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:

  2. jcsd
Know someone interested in this topic? Share a link to this question via email, Google+, Twitter, or Facebook

0
Draft saved Draft deleted