Advanced Search

Question: I'm having problems with plotting graphs of system of equatiosn with control using Gauss Siedel iterative scheme. i want to plot graphs for the costate variables and control function

Question:I'm having problems with plotting graphs of system of equatiosn with control using Gauss Siedel iterative scheme. i want to plot graphs for the costate variables and control function

ogunmiloro 120 Maple
ode linearalgebra numeric differential-equations plotting gauss-seidel-method
June 08 2019
0


 

> 

NULL
> 

restart; lambda := 0.5e6; A__h := 100; `β__1` := 0.34e-1; `ψ__o` := 0.25e-1; mu := 0.4e-3; theta := .7902; alpha := .11; k := 0.136e-3; z := 0.5e-1; delta := .7; eta := .134; xi := .12; beta := 0.1e-1; Q1 := 1; Q2 := 1; h := .1; n := 100; lambda[1][n] := 0; lambda[2][n] := 0; lambda[3][n] := 0; `λ__4`[n] := 0; u[o][0] := 0
> 

for i from 0 to n-1 do B[i+1] := h*(theta*E[i]+A__h+B[i])/(1+`ψ__o`*h*P[i]/(P[i]+xi)+mu+u[o][i]+beta*C[i]); C[i+1] := h*B[i+1][beta*C[i]-psi[o]*P[i]/(P[i]+xi)]-(alpha+mu+k)*C[i+1]; P[i+1] := h*z*C[i+1]/(1+delta+eta); E[i+1] := h*(alpha*C[i+1]+B[i+1]*u[o][i]+E[i])/(1+theta+mu); lambda[1][n-k-1] := -h*[Q1+(beta*C[i]+psi[o]*P[i]/(P[i]+xi))*lambda[2][n-k]-lambda[1][n-k]-mu+u[o][i]]/(1+beta*C[i]+P[i]/(P[i]+xi)); lambda[2][n-k-1] := -h*[-beta*B[i+1]*lambda[1][n-k-1]+Q2-alpha-mu-kappa+z*lambda[3][n-k]+alpha*lambda[4][n-k]-lambda[2][n-k]]/(-beta*B[i+1]+1); lambda[3][n-k-1] := h*[(psi[o]*B[i+1]*P[i+1]/(P[i+1]+xi)^2+theta)*lambda[1][n-k-1]-lambda[2][n-k-1]*psi[o]*B[i+1]*P[i+1]/(P[i+1]+xi)^2-lambda[3][n-k]]/(1+delta+eta); lambda[4][n-k-1] := h*lambda[4][n-k]/(1-theta-mu); R1 := (lambda[1][n-k-1]-lambda[4][n-k-1])*S[i+1]/tau; u[o][i+1] := min(1, max(R1, 0)) end do

Error, recursive assignment

  
> 

seq(i, i = 0 .. 30):
> 

with(plots); with(DEtools); ODE1 := diff(B(T), T) = A__h-beta[1]*B(T)*C(T)-psi[o]*B(T)*P(T)/(P(T)+xi)-mu*B(T)+theta*E(T); ODE2 := diff(C(T), T) = beta[1]*B(T)*C(T)-psi[o]*B(T)*P(T)/(P(T)+xi)-(alpha+mu+k)*C(T); ODE3 := diff(P(T), T) = z*P(T)-(delta+eta)*P(T); ODE4 := diff(E(T), T) = alpha*C(T)-(theta+mu)*E(T); SIRsys := [ODE1, ODE2, ODE3, ODE4]; Bplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, B = 0.167e9 .. 0.17e9, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, B(T)], thickness = 2, linecolor = red, stepsize = .1); Cplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, C = 0 .. 400, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, C(t)], thickness = 2, linecolor = red, stepsize = .1); Pplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, P = 0 .. 0.1e4, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, P(T)], thickness = 2, linecolor = red, stepsize = .1); Eplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, E = 0 .. 0.1e4, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, E(t)], thickness = 2, linecolor = red, stepsize = .1); A := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(B[i], i = 0 .. 100))); p1 := pointplot(A, color = blue, symbol = diamond); display([p1, Bplot])

SIRsys := [diff((table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T), T) = 100-0.1e-1[1]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*C(T)-psi[o]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*P(T)/(P(T)+.12)-0.4e-3*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)+.7902*E(T), diff(C(T), T) = 0.1e-1[1]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*C(T)-psi[o]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*P(T)/(P(T)+.12)-.110536*C(T), diff(P(T), T) = -.784*P(T), diff(E(T), T) = .11*C(T)-.7906*E(T)]

  

A := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  

Error, (in plots:-display) expecting plot structures but received: [p1]

  
> 

G := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(C[i], i = 0 .. 100))); p2 := pointplot(G, color = blue, symbol = soliddiamond); display([p2, Cplot])

G := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  

Error, (in plots:-display) expecting plot structures but received: [p2]

  
> 

``


 

Download cont_plots.mw
 

> 

NULL
> 

restart; lambda := 0.5e6; A__h := 100; `&beta;__1` := 0.34e-1; `&psi;__o` := 0.25e-1; mu := 0.4e-3; theta := .7902; alpha := .11; k := 0.136e-3; z := 0.5e-1; delta := .7; eta := .134; xi := .12; beta := 0.1e-1; Q1 := 1; Q2 := 1; h := .1; n := 100; lambda[1][n] := 0; lambda[2][n] := 0; lambda[3][n] := 0; `&lambda;__4`[n] := 0; u[o][0] := 0
> 

for i from 0 to n-1 do B[i+1] := h*(theta*E[i]+A__h+B[i])/(1+`&psi;__o`*h*P[i]/(P[i]+xi)+mu+u[o][i]+beta*C[i]); C[i+1] := h*B[i+1][beta*C[i]-psi[o]*P[i]/(P[i]+xi)]-(alpha+mu+k)*C[i+1]; P[i+1] := h*z*C[i+1]/(1+delta+eta); E[i+1] := h*(alpha*C[i+1]+B[i+1]*u[o][i]+E[i])/(1+theta+mu); lambda[1][n-k-1] := -h*[Q1+(beta*C[i]+psi[o]*P[i]/(P[i]+xi))*lambda[2][n-k]-lambda[1][n-k]-mu+u[o][i]]/(1+beta*C[i]+P[i]/(P[i]+xi)); lambda[2][n-k-1] := -h*[-beta*B[i+1]*lambda[1][n-k-1]+Q2-alpha-mu-kappa+z*lambda[3][n-k]+alpha*lambda[4][n-k]-lambda[2][n-k]]/(-beta*B[i+1]+1); lambda[3][n-k-1] := h*[(psi[o]*B[i+1]*P[i+1]/(P[i+1]+xi)^2+theta)*lambda[1][n-k-1]-lambda[2][n-k-1]*psi[o]*B[i+1]*P[i+1]/(P[i+1]+xi)^2-lambda[3][n-k]]/(1+delta+eta); lambda[4][n-k-1] := h*lambda[4][n-k]/(1-theta-mu); R1 := (lambda[1][n-k-1]-lambda[4][n-k-1])*S[i+1]/tau; u[o][i+1] := min(1, max(R1, 0)) end do

Error, recursive assignment

  
> 

seq(i, i = 0 .. 30):
> 

with(plots); with(DEtools); ODE1 := diff(B(T), T) = A__h-beta[1]*B(T)*C(T)-psi[o]*B(T)*P(T)/(P(T)+xi)-mu*B(T)+theta*E(T); ODE2 := diff(C(T), T) = beta[1]*B(T)*C(T)-psi[o]*B(T)*P(T)/(P(T)+xi)-(alpha+mu+k)*C(T); ODE3 := diff(P(T), T) = z*P(T)-(delta+eta)*P(T); ODE4 := diff(E(T), T) = alpha*C(T)-(theta+mu)*E(T); SIRsys := [ODE1, ODE2, ODE3, ODE4]; Bplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, B = 0.167e9 .. 0.17e9, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, B(T)], thickness = 2, linecolor = red, stepsize = .1); Cplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, C = 0 .. 400, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, C(t)], thickness = 2, linecolor = red, stepsize = .1); Pplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, P = 0 .. 0.1e4, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, P(T)], thickness = 2, linecolor = red, stepsize = .1); Eplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, E = 0 .. 0.1e4, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, E(t)], thickness = 2, linecolor = red, stepsize = .1); A := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(B[i], i = 0 .. 100))); p1 := pointplot(A, color = blue, symbol = diamond); display([p1, Bplot])

SIRsys := [diff((table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T), T) = 100-0.1e-1[1]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*C(T)-psi[o]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*P(T)/(P(T)+.12)-0.4e-3*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)+.7902*E(T), diff(C(T), T) = 0.1e-1[1]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*C(T)-psi[o]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*P(T)/(P(T)+.12)-.110536*C(T), diff(P(T), T) = -.784*P(T), diff(E(T), T) = .11*C(T)-.7906*E(T)]

  

A := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  

Error, (in plots:-display) expecting plot structures but received: [p1]

  
> 

G := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(C[i], i = 0 .. 100))); p2 := pointplot(G, color = blue, symbol = soliddiamond); display([p2, Cplot])

G := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  

Error, (in plots:-display) expecting plot structures but received: [p2]

  
> 

``


 

Download cont_plots.mw
 

> 

NULL
> 

restart; lambda := 0.5e6; A__h := 100; `&beta;__1` := 0.34e-1; `&psi;__o` := 0.25e-1; mu := 0.4e-3; theta := .7902; alpha := .11; k := 0.136e-3; z := 0.5e-1; delta := .7; eta := .134; xi := .12; beta := 0.1e-1; Q1 := 1; Q2 := 1; h := .1; n := 100; lambda[1][n] := 0; lambda[2][n] := 0; lambda[3][n] := 0; `&lambda;__4`[n] := 0; u[o][0] := 0
> 

for i from 0 to n-1 do B[i+1] := h*(theta*E[i]+A__h+B[i])/(1+`&psi;__o`*h*P[i]/(P[i]+xi)+mu+u[o][i]+beta*C[i]); C[i+1] := h*B[i+1][beta*C[i]-psi[o]*P[i]/(P[i]+xi)]-(alpha+mu+k)*C[i+1]; P[i+1] := h*z*C[i+1]/(1+delta+eta); E[i+1] := h*(alpha*C[i+1]+B[i+1]*u[o][i]+E[i])/(1+theta+mu); lambda[1][n-k-1] := -h*[Q1+(beta*C[i]+psi[o]*P[i]/(P[i]+xi))*lambda[2][n-k]-lambda[1][n-k]-mu+u[o][i]]/(1+beta*C[i]+P[i]/(P[i]+xi)); lambda[2][n-k-1] := -h*[-beta*B[i+1]*lambda[1][n-k-1]+Q2-alpha-mu-kappa+z*lambda[3][n-k]+alpha*lambda[4][n-k]-lambda[2][n-k]]/(-beta*B[i+1]+1); lambda[3][n-k-1] := h*[(psi[o]*B[i+1]*P[i+1]/(P[i+1]+xi)^2+theta)*lambda[1][n-k-1]-lambda[2][n-k-1]*psi[o]*B[i+1]*P[i+1]/(P[i+1]+xi)^2-lambda[3][n-k]]/(1+delta+eta); lambda[4][n-k-1] := h*lambda[4][n-k]/(1-theta-mu); R1 := (lambda[1][n-k-1]-lambda[4][n-k-1])*S[i+1]/tau; u[o][i+1] := min(1, max(R1, 0)) end do

Error, recursive assignment

  
> 

seq(i, i = 0 .. 30):
> 

with(plots); with(DEtools); ODE1 := diff(B(T), T) = A__h-beta[1]*B(T)*C(T)-psi[o]*B(T)*P(T)/(P(T)+xi)-mu*B(T)+theta*E(T); ODE2 := diff(C(T), T) = beta[1]*B(T)*C(T)-psi[o]*B(T)*P(T)/(P(T)+xi)-(alpha+mu+k)*C(T); ODE3 := diff(P(T), T) = z*P(T)-(delta+eta)*P(T); ODE4 := diff(E(T), T) = alpha*C(T)-(theta+mu)*E(T); SIRsys := [ODE1, ODE2, ODE3, ODE4]; Bplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, B = 0.167e9 .. 0.17e9, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, B(T)], thickness = 2, linecolor = red, stepsize = .1); Cplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, C = 0 .. 400, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, C(t)], thickness = 2, linecolor = red, stepsize = .1); Pplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, P = 0 .. 0.1e4, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, P(T)], thickness = 2, linecolor = red, stepsize = .1); Eplot := ODEplot(SIRsys, [B(T), C(T), P(T), E(T)], T = 0 .. 100, E = 0 .. 0.1e4, [[B(0) = 0.17e9, C(0) = 0, P(0) = 400, E(0) = 20]], scene = [T, E(t)], thickness = 2, linecolor = red, stepsize = .1); A := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(B[i], i = 0 .. 100))); p1 := pointplot(A, color = blue, symbol = diamond); display([p1, Bplot])

SIRsys := [diff((table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T), T) = 100-0.1e-1[1]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*C(T)-psi[o]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*P(T)/(P(T)+.12)-0.4e-3*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)+.7902*E(T), diff(C(T), T) = 0.1e-1[1]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*C(T)-psi[o]*(table([1 = (.1*(B[0]+.7902*E[0]+100))/(1.0004+0.25e-2*P[0]/(P[0]+.12)+0.1e-1*C[0])]))(T)*P(T)/(P(T)+.12)-.110536*C(T), diff(P(T), T) = -.784*P(T), diff(E(T), T) = .11*C(T)-.7906*E(T)]

  

A := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  

Error, (in plots:-display) expecting plot structures but received: [p1]

  
> 

G := `<|>`(`<,>`(seq(i, i = 0 .. 100)), `<,>`(seq(C[i], i = 0 .. 100))); p2 := pointplot(G, color = blue, symbol = soliddiamond); display([p2, Cplot])

G := Vector(4, {(1) = ` 101 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

  

Error, (in plots:-pointplot) points cannot be converted to floating-point values

  

Error, (in plots:-display) expecting plot structures but received: [p2]

  
> 

``


 

Download cont_plots.mw

 

 
Please Wait...