Items tagged with numeric numeric Tagged Items Feed

 

Hi

Please find the attachment

It seems a singularity has been occurred at left end point for n less than unit (say n=0.5, n>0)

Is there a way to fix it?

 

n.mw

DEAR SIR,

PLEASE HELP ME WITH THAT QUESTION

 

DEAR SIR

ANYONE CAN HELP TO COMPUTE TIME IN DSOLVE COMMAND?

restart

with(plots)

Nb := 0.1e-4; Nt := 0.1e-4; Sc := 3.0; Sb := 15; Pe := 1; Bi := .5; Pr := 6.8; c[4] := 0; c[6] := .3; c[8] := .4; k[1] := 0; k[2] := 1; k[3] := 0; Un := .1; M := .5

Eq1 := (101-100*lambda)*(1+c[2]*phi(eta))*(diff(f(eta), `$`(eta, 3)))+(diff(f(eta), `$`(eta, 2)))*(f(eta)+g(eta)+c[2]*(diff(phi(eta), eta)))-(diff(f(eta), eta))^2-M*(diff(f(eta), eta)-k[2])+k[2]+Un*(k[2]-(1/2)*eta*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))) = 0

(101-100*lambda)*(1+c[2]*phi(eta))*(diff(diff(diff(f(eta), eta), eta), eta))+(diff(diff(f(eta), eta), eta))*(f(eta)+g(eta)+c[2]*(diff(phi(eta), eta)))-(diff(f(eta), eta))^2-.6*(diff(f(eta), eta))+1.6-0.5000000000e-1*eta*(diff(diff(f(eta), eta), eta)) = 0

(1)

Eq2 := (101-100*lambda)*(1+c[2]*phi(eta))*(diff(g(eta), `$`(eta, 3)))+(diff(g(eta), `$`(eta, 2)))*(f(eta)+g(eta)+c[2]*(diff(phi(eta), eta)))-(diff(g(eta), eta))^2-M*(diff(g(eta), eta)-k[2])+k[2]+Un*(k[2]-(1/2)*eta*(diff(g(eta), `$`(eta, 2)))-(diff(g(eta), eta))) = 0

(101-100*lambda)*(1+c[2]*phi(eta))*(diff(diff(diff(g(eta), eta), eta), eta))+(diff(diff(g(eta), eta), eta))*(f(eta)+g(eta)+c[2]*(diff(phi(eta), eta)))-(diff(g(eta), eta))^2-.6*(diff(g(eta), eta))+1.6-0.5000000000e-1*eta*(diff(diff(g(eta), eta), eta)) = 0

(2)

Eq3 := (1+c[4]*phi(eta))*(diff(theta(eta), `$`(eta, 2)))+Pr*(diff(theta(eta), eta))*(f(eta)+g(eta))+Nb*Pr*(diff(theta(eta), eta))*(diff(phi(eta), eta))*(1+c[6]*(2*phi(eta)+1))+Nt*Pr*(diff(theta(eta), eta))^2+c[4]*(diff(theta(eta), eta))*(diff(phi(eta), eta))-(1/2)*Pr*eta*Un*(diff(theta(eta), eta)) = 0

diff(diff(theta(eta), eta), eta)+6.8*(diff(theta(eta), eta))*(f(eta)+g(eta))+0.68e-4*(diff(theta(eta), eta))*(diff(phi(eta), eta))*(1.3+.6*phi(eta))+0.68e-4*(diff(theta(eta), eta))^2-.3400000000*eta*(diff(theta(eta), eta)) = 0

(3)

Eq4 := (1+c[6]*phi(eta))*(diff(phi(eta), `$`(eta, 2)))+Sc*(f(eta)+g(eta))*(diff(phi(eta), eta))+c[6]*(diff(phi(eta), eta))^2+Nt*(diff(theta(eta), `$`(eta, 2)))/Nb-(1/2)*Sc*eta*Un*(diff(phi(eta), eta)) = 0

(1+.3*phi(eta))*(diff(diff(phi(eta), eta), eta))+3.0*(f(eta)+g(eta))*(diff(phi(eta), eta))+.3*(diff(phi(eta), eta))^2+1.000000000*(diff(diff(theta(eta), eta), eta))-.1500000000*eta*(diff(phi(eta), eta)) = 0

(4)

Eq5 := (1+c[8]*phi(eta))*(diff(chi(eta), `$`(eta, 2)))+c[8]*(diff(phi(eta), eta))*(diff(chi(eta), eta))-Pe*(chi(eta)*(diff(phi(eta), `$`(eta, 2)))+(diff(phi(eta), eta))*(diff(chi(eta), eta)))+Sb*(diff(chi(eta), eta))*(f(eta)+g(eta))-(1/2)*Sb*eta*Un*(diff(chi(eta), eta)) = 0

(1+.4*phi(eta))*(diff(diff(chi(eta), eta), eta))-.6*(diff(phi(eta), eta))*(diff(chi(eta), eta))-chi(eta)*(diff(diff(phi(eta), eta), eta))+15*(diff(chi(eta), eta))*(f(eta)+g(eta))-.7500000000*eta*(diff(chi(eta), eta)) = 0

(5)

Vc[2] := [.2, .3, .4]

etainf := 1.85

bcs := (D(f))(0) = k[1], (D(g))(0) = k[3], f(0) = 0, g(0) = 0, (D(theta))(0) = -Bi*(1-theta(0))/(1+c[4]*phi(0)), Nb*(D(phi))(0)*(1+c[6]*(2*phi(0)+1))+Nt*(D(theta))(0) = 0, chi(0) = 1, (D(f))(etainf) = k[2], (D(g))(etainf) = k[2], theta(etainf) = 0, phi(etainf) = 0, chi(etainf) = 0

(D(f))(0) = 0, (D(g))(0) = 0, f(0) = 0, g(0) = 0, (D(theta))(0) = -.5+.5*theta(0), 0.1e-4*(D(phi))(0)*(1.3+.6*phi(0))+0.1e-4*(D(theta))(0) = 0, chi(0) = 1, (D(f))(1.85) = 1, (D(g))(1.85) = 1, theta(1.85) = 0, phi(1.85) = 0, chi(1.85) = 0

(6)

dsys := {Eq1, Eq2, Eq3, Eq4, Eq5, bcs}

for i to 3 do c[2] := Vc[2][i]; dsol[i] := dsolve(dsys, numeric, continuation = lambda); print(c[2]); print(dsol[i](0)) end do

.2

 

[eta = 0., chi(eta) = HFloat(1.0), diff(chi(eta), eta) = HFloat(3.888290578689045), f(eta) = HFloat(0.0), diff(f(eta), eta) = HFloat(0.0), diff(diff(f(eta), eta), eta) = HFloat(2.244199282192492), g(eta) = HFloat(0.0), diff(g(eta), eta) = HFloat(0.0), diff(diff(g(eta), eta), eta) = HFloat(2.244199282192492), phi(eta) = HFloat(-2.044191234673432), diff(phi(eta), eta) = HFloat(5.227515304629519), theta(eta) = HFloat(0.2317093657771352), diff(theta(eta), eta) = HFloat(-0.38414531711143246)]

 

.3

 

[eta = 0., chi(eta) = HFloat(1.0), diff(chi(eta), eta) = HFloat(4.148187853914835), f(eta) = HFloat(0.0), diff(f(eta), eta) = HFloat(0.0), diff(diff(f(eta), eta), eta) = HFloat(3.1086884419918364), g(eta) = HFloat(0.0), diff(g(eta), eta) = HFloat(0.0), diff(diff(g(eta), eta), eta) = HFloat(3.1086884419918364), phi(eta) = HFloat(-2.049332060722701), diff(phi(eta), eta) = HFloat(5.527786294980874), theta(eta) = HFloat(0.2216792480031605), diff(theta(eta), eta) = HFloat(-0.38916037599841974)]

 

.4

 

[eta = 0., chi(eta) = HFloat(0.9999999999999998), diff(chi(eta), eta) = HFloat(4.580796631072469), f(eta) = HFloat(0.0), diff(f(eta), eta) = HFloat(0.0), diff(diff(f(eta), eta), eta) = HFloat(5.687607599246298), g(eta) = HFloat(0.0), diff(g(eta), eta) = HFloat(0.0), diff(diff(g(eta), eta), eta) = HFloat(5.687607599246298), phi(eta) = HFloat(-2.0568809171520708), diff(phi(eta), eta) = HFloat(6.0203396482123575), theta(eta) = HFloat(0.20686299926628238), diff(theta(eta), eta) = HFloat(-0.39656850036685887)]

(7)

NULL

 

 

Download 3DAKc2w_-_Copy.mw

hi .please help me for solve this equations.

bbb2.mw

restart; d[11] := 1; mu[11] := 1; q[311] := 1; d[33] := 1; mu[33] := 1; a[11] := 1; e[311] := 1; a[33] := 1; A := 1; g[111111] := 1; c[1111] := 1; g[113113] := 1; f[3113] := 1; beta[11] := 1; `ΔT` := 1; II := 1; L := 1

J := d[11]*(diff(Phi(x, z), x, x))+mu[11]*(diff(psi(x, z), x, x))+q[311]*(diff(w(x), x, x))+d[33]*(diff(Phi(x, z), z, z))+mu[33]*(diff(psi(x, z), z, z));

diff(diff(Phi(x, z), x), x)+diff(diff(psi(x, z), x), x)+diff(diff(w(x), x), x)+diff(diff(Phi(x, z), z), z)+diff(diff(psi(x, z), z), z)

(1)

B := a[11]*(diff(Phi(x, z), x, x))+d[11]*(diff(psi(x, z), x, x))+e[311]*(diff(w(x), x, x))+a[33]*(diff(Phi(x, z), z, z))+d[33]*(diff(psi(x, z), z, z));

diff(diff(Phi(x, z), x), x)+diff(diff(psi(x, z), x), x)+diff(diff(w(x), x), x)+diff(diff(Phi(x, z), z), z)+diff(diff(psi(x, z), z), z)

(2)

R := A*(g[111111]*(diff(u[0](x), x, x, x, x))-c[1111]*(diff(u[0](x), x, x)+(1/2)*(diff((diff(w(x), x))^2, x)))+e[311]*(diff(diff(Phi(x, z), z), x))+q[311]*(diff(diff(psi(x, z), z), x)));

diff(diff(diff(diff(u[0](x), x), x), x), x)-(diff(diff(u[0](x), x), x))-(diff(w(x), x))*(diff(diff(w(x), x), x))+diff(diff(Phi(x, z), x), z)+diff(diff(psi(x, z), x), z)

(3)

S := -II*g[111111]*(diff(w(x), x, x, x, x, x, x))-II*c[1111]*(diff(w(x), x, x, x, x))+A*g[113113]*(diff(w(x), x, x, x, x))-A*f[3113]*(diff(diff(Phi(x, z), z), x, x))-A*(c[1111]*(diff(u[0](x), x, x)+(1/2)*(diff((diff(w(x), x))^2, x)))+e[311]*(diff(diff(Phi(x, z), z), x))+q[311]*(diff(diff(psi(x, z), z), x)))*(diff(w(x), x))-A*(diff(w(x), x, x))*(c[1111]*(diff(u[0](x), x)+(1/2)*(diff(w(x), x))^2)+e[311]*(diff(Phi(x, z), z))+q[311]*(diff(psi(x, z), z))-beta[11]*`ΔT`);

-(diff(diff(diff(diff(diff(diff(w(x), x), x), x), x), x), x))-(diff(diff(diff(Phi(x, z), x), x), z))-(diff(diff(u[0](x), x), x)+(diff(w(x), x))*(diff(diff(w(x), x), x))+diff(diff(Phi(x, z), x), z)+diff(diff(psi(x, z), x), z))*(diff(w(x), x))-(diff(diff(w(x), x), x))*(diff(u[0](x), x)+(1/2)*(diff(w(x), x))^2+diff(Phi(x, z), z)+diff(psi(x, z), z)-1)

(4)

dsys := {B, J, R, S}; BCS := {D@@2*w(0) = 0, D@@2*w(L) = 0, Phi(x = 0) = 0, Phi(x = L) = 0, Phi(z = -(1/2)*h) = 0, Phi(z = (1/2)*h) = 0, psi(x = 0) = 0, psi(x = L) = 0, psi(z = -(1/2)*h) = 0, psi(z = (1/2)*h) = 0, w(x = 0) = 0, w(x = L) = 0, u[0](x = 0) = 0, u[0](x = L) = 0, (D(w))(0) = 0, (D(w))(L) = 0, (D(u[0]))(0) = 0, (D(u[0]))(L) = 0}

{D@@2*w(0) = 0, D@@2*w(L) = 0, Phi(x = 0) = 0, Phi(x = L) = 0, Phi(z = -(1/2)*h) = 0, Phi(z = (1/2)*h) = 0, psi(x = 0) = 0, psi(x = L) = 0, psi(z = -(1/2)*h) = 0, psi(z = (1/2)*h) = 0, w(x = 0) = 0, w(x = L) = 0, u[0](x = 0) = 0, u[0](x = L) = 0, (D(w))(0) = 0, (D(w))(L) = 0, (D(u[0]))(0) = 0, (D(u[0]))(L) = 0}

(5)

dsol5 := dsolve(dsys, numeric)

Error, (in dsolve/numeric/process_input) missing differential equations and initial or boundary conditions in the first argument: dsys

 

NULL

NULL

NULL

if former equations are not solvable , please help me for another way, in which at first two equation solve..in this way in equation [J and B] assume that q[311]=e[311]=0 and dsolve perform to find Φ and  ψ

after by finding Φ and  ψ is use for detemine w and u0

please see attached file below[bbb2_2.mw]

bbb2_2.mw

Download bbb2.mw

Dear All,

I am going to solve the following systems of ODEs but get the error: Newton iteration is not converging.
Could you please share your idea with me. In the case of AA=-0.2,0,0.2,0.4,...; I could get the solution.
Thank you in advance.


restart;
with(plots);
Pr := 2; Le := 2; nn := 2; Nb := .1; Nt := .1; QQ := .1; SS := .1; BB := .1; CC := .1; Ec := .1; MM := .2;AA:=-0.4;

Eq1 := diff(f(eta), `$`(eta, 3))+f(eta).(diff(f(eta), `$`(eta, 2)))-2.*nn/(nn+1).((diff(f(eta), eta))^2)-MM.(diff(f(eta), eta)) = 0; Eq2 := 1/Pr.(diff(theta(eta), `$`(eta, 2)))+f(eta).(diff(theta(eta), eta))-4.*nn/(nn+1).(diff(f(eta), eta)).theta(eta)+Nb.(diff(theta(eta), eta)).(diff(h(eta), eta))+Nt.((diff(theta(eta), eta))^2)+Ec.((diff(f(eta), `$`(eta, 2)))^2)-QQ.theta(eta) = 0;
Eq3 := diff(h(eta), `$`(eta, 2))+Le.f(eta).(diff(h(eta), eta))+Nt/Nb.(diff(theta(eta), `$`(eta, 2))) = 0;

bcs := f(0) = SS, (D(f))(0) = 1+AA.((D@@2)(f))(0), theta(0) = 1+BB.(D(theta))(0), phi(0) = 1+CC.(D(phi))(0), (D(f))(etainf) = 0, theta(etainf) = 0, phi(etainf) = 0

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

I've been trying to numerically solve (and plot) this equation. Maple tells me that some matrix is singular - I have no idea, what I can do.

eq := diff(y(x), `$`(x, 3))-(diff(y(x), x))*(diff(y(x), x))+1 = 0;

cond := (D(y))(0) = 0, (D(y))(1) = 1, ((D@@2)(y))(0) = 0

de := dsolve({cond, eq}, y(x), numeric);

Error, (in dsolve/numeric/bvp) matrix is singular

I am unable to solve the attached optimal control problem,please any one who many help  me in guideing .tnx

restart:
unprotect('gamma');

L:=b[1]*c(t)+b[2]*i(t)+w[1]*(u[1])^2/2+w[2]*(u[2])^2/2+w[3]*(u[3])^2/2;
1 2 1 2 1 2
b[1] c(t) + b[2] i(t) + - w[1] u[1] + - w[2] u[2] + - w[3] u[3]
2 2 2
H:=L+lambda[1](t)*((1-p*Psi)*tau+phi* v + delta *r-lambda*(1-u[3])*s-u[1]*varphi*s -mu*s ) +lambda[2](t)*(p*Psi*tau + u[1]*vartheta*s -gamma*lambda* (1-u[3])*v-(mu+phi)*v ) +lambda[3](t)*( (1-u[3])*rho*lambda* (s +gamma*v)+(1-q)* u[2]*eta*i -(mu +beta +chi)*c ) +lambda[4](t)* ((1-rho)*(1-u[3])*lambda*( s +gamma*v) +chi*c - u[2]*eta*i - (mu +alpha )*i) +lambda[5](t)*( beta*c + u[2]*q*eta*i -(mu +delta)*r);
1 2 1 2 1 2
b[1] c(t) + b[2] i(t) + - w[1] u[1] + - w[2] u[2] + - w[3] u[3] + lambda[1](t
2 2 2

) ((1 - p Psi) tau + phi v + delta r - lambda (1 - u[3]) s - u[1] varphi s

- mu s) + lambda[2](t) (p Psi tau + u[1] vartheta s

- gamma lambda (1 - u[3]) v - (mu + phi) v) + lambda[3](t) ((1 - u[3]) rho

lambda (s + gamma v) + (1 - q) u[2] eta i - (mu + beta + chi) c) + lambda[4](t

) ((1 - rho) (1 - u[3]) lambda (s + gamma v) + chi c - u[2] eta i

- (mu + alpha) i) + lambda[5](t) (beta c + u[2] q eta i - (mu + delta) r)
du1:=diff(H,u[1]);

w[1] u[1] - lambda[1](t) varphi s + lambda[2](t) vartheta s
du2:=diff(H,u[2]);du3:=diff(H,u[3]);
w[2] u[2] + lambda[3](t) (1 - q) eta i - lambda[4](t) eta i

+ lambda[5](t) q eta i
w[3] u[3] + lambda[1](t) lambda s + lambda[2](t) gamma lambda v

- lambda[3](t) rho lambda (s + gamma v)

- lambda[4](t) (1 - rho) lambda (s + gamma v)

ddu1 := -A[1] u[1] + psi[1](t) beta x[1] x[3] - psi[2](t) beta x[1] x[3]

ddu2 := -A[2] u[2] - psi[3](t) k x[2]
sol_u1 := solve(du1, u[1]);
s(t) (lambda[1](t) varphi - lambda[2](t) vartheta)
--------------------------------------------------
w[1]
sol_u2 := solve(du2, u[2]);sol_u3 := solve(du3, u[3]);
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q)
----------------------------------------------------------------------
w[2]
1
---- (lambda (-lambda[1](t) s - lambda[2](t) gamma v + lambda[3](t) rho s
w[3]

+ lambda[3](t) rho gamma v + lambda[4](t) s + lambda[4](t) gamma v

- lambda[4](t) rho s - lambda[4](t) rho gamma v))
Dx2:=subs(u[1]= s*(lambda[1](t)*varphi-lambda[2](t)*vartheta)/w[1] ,u[2]= eta*i*(-lambda[3](t)+lambda[3](t)*q+lambda[4](t)-lambda[5](t)*q)/w[2], u[3]=-lambda*(lambda[1](t)*s+lambda[2](t)*gamma*v-lambda[3](t)*rho*s-lambda[3](t)*rho*gamma*v-lambda[4](t)*s-lambda[4](t)*gamma*v+lambda[4](t)*rho*s+lambda[4](t)*rho*gamma*v)/w[3] ,H );
2 2
s (lambda[1](t) varphi - lambda[2](t) vartheta)
b[1] c(t) + b[2] i(t) + -------------------------------------------------
2 w[1]

2 2 2
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q)
+ ------------------------------------------------------------------------- +
2 w[2]

1 / 2
------ \lambda (lambda[1](t) s + lambda[2](t) gamma v - lambda[3](t) rho s
2 w[3]

- lambda[3](t) rho gamma v - lambda[4](t) s - lambda[4](t) gamma v

/
\ |
+ lambda[4](t) rho s + lambda[4](t) rho gamma v)^2/ + lambda[1](t) |(1
\

/ 1
- p Psi) tau + phi v + delta r - lambda |1 + ---- (lambda (lambda[1](t) s
\ w[3]

+ lambda[2](t) gamma v - lambda[3](t) rho s - lambda[3](t) rho gamma v

- lambda[4](t) s - lambda[4](t) gamma v + lambda[4](t) rho s

\
+ lambda[4](t) rho gamma v))| s
/

2 \
s (lambda[1](t) varphi - lambda[2](t) vartheta) varphi |
- ------------------------------------------------------- - mu s| +
w[1] /

/
|
lambda[2](t) |p Psi tau
\

2
s (lambda[1](t) varphi - lambda[2](t) vartheta) vartheta /
+ --------------------------------------------------------- - gamma lambda |1 +
w[1] \

1
---- (lambda (lambda[1](t) s + lambda[2](t) gamma v - lambda[3](t) rho s
w[3]

- lambda[3](t) rho gamma v - lambda[4](t) s - lambda[4](t) gamma v

\
\ |
+ lambda[4](t) rho s + lambda[4](t) rho gamma v))| v - (mu + phi) v| +
/ /

// 1
lambda[3](t) ||1 + ---- (lambda (lambda[1](t) s + lambda[2](t) gamma v
\\ w[3]

- lambda[3](t) rho s - lambda[3](t) rho gamma v - lambda[4](t) s

\
- lambda[4](t) gamma v + lambda[4](t) rho s + lambda[4](t) rho gamma v))|
/

1 / 2 2
rho lambda (s + gamma v) + ---- \(1 - q) eta i (-lambda[3](t)
w[2]

\ \
+ lambda[3](t) q + lambda[4](t) - lambda[5](t) q)/ - (mu + beta + chi) c| +
/

/
| / 1
lambda[4](t) |(1 - rho) |1 + ---- (lambda (lambda[1](t) s
\ \ w[3]

+ lambda[2](t) gamma v - lambda[3](t) rho s - lambda[3](t) rho gamma v

- lambda[4](t) s - lambda[4](t) gamma v + lambda[4](t) rho s

\
+ lambda[4](t) rho gamma v))| lambda (s + gamma v) + chi c
/

2 2
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q)
- ------------------------------------------------------------------------
w[2]

\ /
| |
- (mu + alpha) i| + lambda[5](t) |beta c
/ \

+

2 2
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q) q
--------------------------------------------------------------------------
w[2]

\
|
- (mu + delta) r|
/
ode1:=diff(lambda[1](t),t)=-diff(H,s);ode2:=diff(lambda[2](t),t)=-diff(H,v);ode3:=diff(psi[3](t),t)=-diff(H,c);ode4:=diff(lambda[4](t),t)=-diff(H,i);ode5:=diff(lambda[5](t),t)=-diff(H,r);
d
--- lambda[1](t) = -lambda[1](t) (-lambda (1 - u[3]) - u[1] varphi - mu)
dt

- lambda[2](t) u[1] vartheta - lambda[3](t) (1 - u[3]) rho lambda

- lambda[4](t) (1 - rho) (1 - u[3]) lambda
d
--- lambda[2](t) = -lambda[1](t) phi
dt

- lambda[2](t) (-gamma lambda (1 - u[3]) - mu - phi)

- lambda[3](t) (1 - u[3]) rho lambda gamma

- lambda[4](t) (1 - rho) (1 - u[3]) lambda gamma
d
--- psi[3](t) = -lambda[3](t) (-mu - beta - chi) - lambda[4](t) chi
dt

- lambda[5](t) beta
d
--- lambda[4](t) = -lambda[3](t) (1 - q) u[2] eta
dt

- lambda[4](t) (-u[2] eta - mu - alpha) - lambda[5](t) u[2] q eta
d
--- lambda[5](t) = -lambda[1](t) delta - lambda[5](t) (-mu - delta)
dt
restart:
#Digits:=10:


unprotect('gamma');
lambda:=0.51:
mu:=0.002:
beta:=0.0115:
delta:=0.003:
alpha:=0.33:
chi:=0.00274:
k:=6.24:
gamma:=0.4:
rho:=0.338:;tau=1000:;Psi:=0.1:;p:=0.6:;phi:=0.001:;eta:=0.001124:q:=0.6:varphi:=0.9:;vatheta:=0.9:
b[1]:=2:;b[2]:=3:;w[1]:=4:;w[2]:=5:;w[3]:=6:
#u[1]:=s(t)*(lambda[1](t)*varphi-lambda[2](t)*vartheta)/w[1]:
#u[2]:=eta*i*(-lambda[3](t)+lambda[3](t)*q+lambda[4](t)-lambda[5](t)*q)/w[2]:;u[3]:=lambda*(-lambda[1](t)*s-lambda[2](t)*gamma*v+lambda[3](t)*rho*s+lambda[3](t)*rho*gamma*v+lambda[4](t)*s+lambda[4](t)*gamma*v-lambda[4](t)*rho*s-lambda[4](t)*rho*gamma*v)/w[3]:
ics := s(0)=8200, v(0)=2800,c(0)=1100,i(0)=1500,r(0)=200,lambda[1](20)=0,lambda[2](20)=0,lambda[3](20)=0,lambda[4](20)=0,lambda[5](20)=0:
ode1:=diff(s(t),t)=(1-p*Psi)*tau+phi* v(t) + delta *r(t)-lambda*(1-u[3])*s(t)-u[1]*varphi*s(t) -mu*s(t),
diff(v(t), t) =p*Psi*tau + u[1]*vartheta*s(t) -gamma*lambda* (1-u[3])*v(t)-(mu+phi)*v(t) ,
diff(c(t), t) =(1-u[3])*rho*lambda* (s(t) +gamma*v(t))+(1-q)* u[2]*eta*i(t) -(mu +beta +chi)*c(t),
diff(i(t), t) =(1-rho)*(1-u[3])*lambda*( s(t) +gamma*v(t)) +chi*c(t) - u[2]*eta*i(t) - (mu +alpha )*i(t),
diff(r(t), t) = beta*c(t) + u[2]*q*eta*i(t) -(mu +delta)*r(t),
diff(lambda[1](t), t) = -lambda[1](t)*(-lambda*(1-u[3])-u[1]*varphi-mu)-lambda[2](t)*u[1]*vartheta-lambda[3](t)*(1-u[3])*rho*lambda-lambda[4](t)*(1-rho)*(1-u[3])*lambda,diff(lambda[2](t),t)=-lambda[1](t)*phi-lambda[2](t)*(-gamma*lambda*(1-u[3])-mu-phi)-lambda[3](t)*(1-u[3])*rho*lambda*gamma-lambda[4](t)*(1-rho)*(1-u[3])*lambda*gamma,diff(lambda[3](t),t)=-lambda[3](t)*(-mu-beta-chi)-lambda[4](t)*chi-lambda[5](t)*beta,diff(lambda[4](t),t)=-lambda[3](t)*(1-q)*u[2]*eta-lambda[4](t)*(-u[2]*eta-mu-alpha)-lambda[5](t)*u[2]*q*eta,diff(lambda[5](t),t)=-lambda[1](t)*delta-lambda[5](t)*(-mu-delta);
d
--- s(t) = (1 - p Psi) tau + phi v(t) + delta r(t) - lambda (1 - u[3]) s(t)
dt

d
- u[1] varphi s(t) - mu s(t), --- v(t) = p Psi tau + u[1] vartheta s(t)
dt

d
- gamma lambda (1 - u[3]) v(t) - (mu + phi) v(t), --- c(t) = (1 - u[3]) rho lambda
dt

(s(t) + gamma v(t)) + (1 - q) u[2] eta - (mu + beta + chi) c(t), 0 = (1

- rho) (1 - u[3]) lambda (s(t) + gamma v(t)) + chi c(t) - u[2] eta - mu

d d
- alpha, --- r(t) = beta c(t) + u[2] q eta - (mu + delta) r(t), ---
dt dt

lambda[1](t) = -lambda[1](t) (-lambda (1 - u[3]) - u[1] varphi - mu)

- lambda[2](t) u[1] vartheta - lambda[3](t) (1 - u[3]) rho lambda

d
- lambda[4](t) (1 - rho) (1 - u[3]) lambda, --- lambda[2](t) =
dt
-lambda[1](t) phi - lambda[2](t) (-gamma lambda (1 - u[3]) - mu - phi)

- lambda[3](t) (1 - u[3]) rho lambda gamma

d
- lambda[4](t) (1 - rho) (1 - u[3]) lambda gamma, --- lambda[3](t) =
dt
d
-lambda[3](t) (-mu - beta - chi) - lambda[4](t) chi - lambda[5](t) beta, ---
dt

lambda[4](t) = -lambda[3](t) (1 - q) u[2] eta

- lambda[4](t) (-u[2] eta - mu - alpha) - lambda[5](t) u[2] q eta,

d
--- lambda[5](t) = -lambda[1](t) delta - lambda[5](t) (-mu - delta)
dt

sol := dsolve({c(0) = 0, i(0) = 0, r(0) = .1, s(0) = 0, v(0) = 0, diff(c(t), t) = (1-u[3])*rho*lambda*(s(t)+gamma*v(t))+(1-q)*u[2]*eta*i(t)-(mu+beta+chi)*c(t), diff(i(t), t) = (1-rho)*(1-u[3])*lambda*(s(t)+gamma*v(t))+chi*c(t)-u[2]*eta*i(t)-(mu+alpha)*i(t), diff(r(t), t) = beta*c(t)+u[2]*q*eta*i(t)-(mu+delta)*r(t), diff(s(t), t) = (1-p*Psi)*tau+phi*v(t)+delta*r(t)-lambda*(1-u[3])*s(t)-u[1]*varphi*s(t)-mu*s(t), diff(v(t), t) = p*Psi*tau+u[1]*vartheta*s(t)-gamma*lambda*(1-u[3])*v(t)-(mu+phi)*v(t), diff(lambda[1](t), t) = -lambda[1](t)*(-lambda*(1-u[3])-u[1]*varphi-mu)-lambda[2](t)*u[1]*vartheta-lambda[3](t)*(1-u[3])*rho*lambda-lambda[4](t)*(1-rho)*(1-u[3])*lambda, diff(lambda[2](t), t) = -lambda[1](t)*phi-lambda[2](t)*(-gamma*lambda*(1-u[3])-mu-phi)-lambda[3](t)*(1-u[3])*rho*lambda*gamma-lambda[4](t)*(1-rho)*(1-u[3])*lambda*gamma, diff(lambda[3](t), t) = -lambda[3](t)*(-mu-beta-chi)-lambda[4](t)*chi-lambda[5](t)*beta, diff(lambda[4](t), t) = -lambda[3](t)*(1-q)*u[2]*eta-lambda[4](t)*(-u[2]*eta-mu-alpha)-lambda[5](t)*u[2]*q*eta, diff(lambda[5](t), t) = -lambda[1](t)*delta-lambda[5](t)*(-mu-delta), lambda[1](20) = 0, lambda[2](20) = 0, lambda[3](20) = 0, lambda[4](20) = 0, lambda[5](20) = 0}, type = numeric);
Error, (in dsolve/numeric/process_input) invalid specification of initial conditions, got 1 = 0

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

Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations

dsolve[':-interactive']({});
Error, `:=` unexpected
sol:=dsolve([ode1,ics],numeric, method = bvp[midrich],maxmesh=500);
Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations

eq1:=diff(s(t), t)=(1-p*Psi)*tau+phi* v(t) + delta *r(t)-lambda*(1-u[3])*s(t)-u[1]*varphi*s(t) -mu*s(t);
eq2:diff(v(t), t) =p*Psi*tau + u[1]*vartheta*s(t) -gamma*lambda* (1-u[3])*v(t)-(mu+phi)*v(t);
eq3:=diff(c(t), t) =(1-u[3])*rho*lambda* (s(t) +gamma*v(t))+(1-q)* u[2]*eta*i(t) -(mu +beta +chi)*c(t);
eq4:=diff(i(t), t) =(1-rho)*(1-u[3])*lambda*( s(t) +gamma*v(t)) +chi*c(t) - u[2]*eta*i(t) - (mu +alpha )*i(t);
eq5:=diff(r(t), t) = beta*c(t) + u[2]*q*eta*i(t) -(mu +delta)*r(t);

d
--- s(t) = (1 - p Psi) tau + phi v(t) + delta r(t) - lambda (1 - u[3]) s(t)
dt

- u[1] varphi s(t) - mu s(t)
d
--- v(t) = p Psi tau + u[1] vartheta s(t) - gamma lambda (1 - u[3]) v(t)
dt

- (mu + phi) v(t)
d
--- c(t) = (1 - u[3]) rho lambda (s(t) + gamma v(t)) + (1 - q) u[2] eta i(t)
dt

- (mu + beta + chi) c(t)
d
--- i(t) = (1 - rho) (1 - u[3]) lambda (s(t) + gamma v(t)) + chi c(t)
dt

- u[2] eta i(t) - (mu + alpha) i(t)
d
--- r(t) = beta c(t) + u[2] q eta i(t) - (mu + delta) r(t)
dt
eq6:=diff(Q(t),t)=b[1]*c(t)+b[2]*i(t)+w[1]*(u[1])^2/2+w[2]*(u[2])^2/2+w[3]*(u[3])^2/2;
d 1 2 1 2 1 2
--- Q(t) = b[1] c(t) + b[2] i(t) + - w[1] u[1] + - w[2] u[2] + - w[3] u[3]
dt 2 2 2
ics:=s(0)=8200, v(0)=2800,c(0)=1100,i(0)=1500,r(0)=200,Q(0)=6700;
s(0) = 8200, v(0) = 2800, c(0) = 1100, i(0) = 1500, r(0) = 200, Q(0) = 6700
sol0:=dsolve({eq1,eq2,eq3,eq4,eq5,eq6,ics},type=numeric,stiff=true,'parameters'=[u[1],u[2],u[3]],abserr=1e-15,relerr=1e-12,maxfun=0,range=0..50):
Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations
with(plots):
Q0:=6700;
6700
obj:=proc(u)
global sol0,Q0;
local ob1;
try
sol0('parameters'=[u[1],u[2],u[3]]):
ob1:=subs(sol0(20.),Q(t)):
catch :
ob1:=0;
end try;
#ob1:=subs(sol0(20.),Q(t));
if ob1>Q0 then Q0:=ob1;print(Q0,u);end;
ob1;
end proc;
proc(u) ... end;
obj([1,1,1]);
0
obj([3,2.5],4);
0
u0:=Vector(3,[0.,0.,0.],datatype=float[8]);
Vector[column](%id = 85973880)

Q0:=0;
Q0 := 0
with(Optimization);
[ImportMPS, Interactive, LPSolve, LSSolve, Maximize, Minimize, NLPSolve,

QPSolve]
sol2:=NLPSolve(3,obj,initialpoint=u0,method=nonlinearsimplex,maximize,evaluationlimit=100):
sol0('parameters'=[3.18125786060723, 2.36800986932868]);
sol0(parameters = [3.18125786060723, 2.36800986932868])
for i from 1 to 3 do odeplot(sol0,[t,x[i](t)],0..20,thickness=3,axes=boxed);od;
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

 

hello , 

how i can exract value from pdsolve ,i need to use dU(x,R)/dR 

thank you 

 

restart; with(plots)

n := 1/3;

1/3

(1)

Uu := (3*n+1)*(1-R^((n+1)/n))/(n+1);

-(3/2)*R^4+3/2

(2)

eq := Uu*(diff(theta(x, R), x))-4*(diff(R*(diff(theta(x, R), R)), R))/R;

(-(3/2)*R^4+3/2)*(diff(theta(x, R), x))-4*(diff(theta(x, R), R)+R*(diff(diff(theta(x, R), R), R)))/R

(3)

IBC := {theta(0, R) = 1, theta(x, 1) = 0, (D[2](theta))(x, 0) = 0};

{theta(0, R) = 1, theta(x, 1) = 0, (D[2](theta))(x, 0) = 0}

(4)

pds := pdsolve(eq, IBC, numeric);

module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

(5)

U := subs(pds:-value(output = listprocedure), theta(x, R));

proc () local tv, xv, solnproc, stype, ndsol, vals; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; Digits := trunc(evalhf(Digits)); solnproc := proc (tv, xv) local INFO, errest, nd, dvars, dary, daryt, daryx, vals, msg, i, j; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; table( [( "soln_procedures" ) = array( 1 .. 1, [( 1 ) = (18446744074366926358)  ] ) ] ) INFO := table( [( "timestep" ) = 0.500000000000000e-1, ( "IBC" ) = b, ( "spaceidx" ) = 2, ( "fdepvars" ) = [theta(x, R)], ( "dependson" ) = [{1}], ( "eqnords" ) = [[1, 2]], ( "intspace" ) = Matrix(21, 1, {(1, 1) = .0, (2, 1) = .0, (3, 1) = .0, (4, 1) = .0, (5, 1) = .0, (6, 1) = .0, (7, 1) = .0, (8, 1) = .0, (9, 1) = .0, (10, 1) = .0, (11, 1) = .0, (12, 1) = .0, (13, 1) = .0, (14, 1) = .0, (15, 1) = .0, (16, 1) = .0, (17, 1) = .0, (18, 1) = .0, (19, 1) = .0, (20, 1) = .0, (21, 1) = .0}, datatype = float[8], order = C_order), ( "solvec2" ) = Vector(21, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0, (21) = .0}, datatype = float[8]), ( "allocspace" ) = 21, ( "solmat_ne" ) = 0, ( "depords" ) = [[1, 2]], ( "BCS", 1 ) = {[[1, 0, 1], b[1, 0, 1]], [[1, 1, 0], b[1, 1, 0]]}, ( "spacepts" ) = 21, ( "solvec3" ) = Vector(21, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0, (21) = .0}, datatype = float[8]), ( "autonomous" ) = true, ( "vectorproc" ) = proc (v, vp, vpp, t, x, k, h, n, vec) local _s1, _s2, _s3, _s4, _s5, _s6, xi; _s3 := 4*k; _s4 := -3*h^2; _s5 := 2*h*k; _s6 := 2*k*h^2; vec[1] := 0; vec[n] := 0; for xi from 2 to n-1 do _s1 := -vp[xi-1]+vp[xi+1]; _s2 := vp[xi-1]-2*vp[xi]+vp[xi+1]; vec[xi] := (_s4*vp[xi]*x[xi]^5+_s2*_s3*x[xi]-_s4*vp[xi]*x[xi]+_s1*_s5)/(_s6*x[xi]) end do end proc, ( "timeidx" ) = 1, ( "extrabcs" ) = [0], ( "pts", R ) = [0, 1], ( "solvec5" ) = 0, ( "timevar" ) = x, ( "t0" ) = 0, ( "solmat_v" ) = Vector(147, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0, (21) = .0, (22) = .0, (23) = .0, (24) = .0, (25) = .0, (26) = .0, (27) = .0, (28) = .0, (29) = .0, (30) = .0, (31) = .0, (32) = .0, (33) = .0, (34) = .0, (35) = .0, (36) = .0, (37) = .0, (38) = .0, (39) = .0, (40) = .0, (41) = .0, (42) = .0, (43) = .0, (44) = .0, (45) = .0, (46) = .0, (47) = .0, (48) = .0, (49) = .0, (50) = .0, (51) = .0, (52) = .0, (53) = .0, (54) = .0, (55) = .0, (56) = .0, (57) = .0, (58) = .0, (59) = .0, (60) = .0, (61) = .0, (62) = .0, (63) = .0, (64) = .0, (65) = .0, (66) = .0, (67) = .0, (68) = .0, (69) = .0, (70) = .0, (71) = .0, (72) = .0, (73) = .0, (74) = .0, (75) = .0, (76) = .0, (77) = .0, (78) = .0, (79) = .0, (80) = .0, (81) = .0, (82) = .0, (83) = .0, (84) = .0, (85) = .0, (86) = .0, (87) = .0, (88) = .0, (89) = .0, (90) = .0, (91) = .0, (92) = .0, (93) = .0, (94) = .0, (95) = .0, (96) = .0, (97) = .0, (98) = .0, (99) = .0, (100) = .0, (101) = .0, (102) = .0, (103) = .0, (104) = .0, (105) = .0, (106) = .0, (107) = .0, (108) = .0, (109) = .0, (110) = .0, (111) = .0, (112) = .0, (113) = .0, (114) = .0, (115) = .0, (116) = .0, (117) = .0, (118) = .0, (119) = .0, (120) = .0, (121) = .0, (122) = .0, (123) = .0, (124) = .0, (125) = .0, (126) = .0, (127) = .0, (128) = .0, (129) = .0, (130) = .0, (131) = .0, (132) = .0, (133) = .0, (134) = .0, (135) = .0, (136) = .0, (137) = .0, (138) = .0, (139) = .0, (140) = .0, (141) = .0, (142) = .0, (143) = .0, (144) = .0, (145) = .0, (146) = .0, (147) = .0}, datatype = float[8], order = C_order, attributes = [source_rtable = (Matrix(21, 7, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (13, 1) = .0, (13, 2) = .0, (13, 3) = .0, (13, 4) = .0, (13, 5) = .0, (13, 6) = .0, (13, 7) = .0, (14, 1) = .0, (14, 2) = .0, (14, 3) = .0, (14, 4) = .0, (14, 5) = .0, (14, 6) = .0, (14, 7) = .0, (15, 1) = .0, (15, 2) = .0, (15, 3) = .0, (15, 4) = .0, (15, 5) = .0, (15, 6) = .0, (15, 7) = .0, (16, 1) = .0, (16, 2) = .0, (16, 3) = .0, (16, 4) = .0, (16, 5) = .0, (16, 6) = .0, (16, 7) = .0, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = .0, (17, 5) = .0, (17, 6) = .0, (17, 7) = .0, (18, 1) = .0, (18, 2) = .0, (18, 3) = .0, (18, 4) = .0, (18, 5) = .0, (18, 6) = .0, (18, 7) = .0, (19, 1) = .0, (19, 2) = .0, (19, 3) = .0, (19, 4) = .0, (19, 5) = .0, (19, 6) = .0, (19, 7) = .0, (20, 1) = .0, (20, 2) = .0, (20, 3) = .0, (20, 4) = .0, (20, 5) = .0, (20, 6) = .0, (20, 7) = .0, (21, 1) = .0, (21, 2) = .0, (21, 3) = .0, (21, 4) = .0, (21, 5) = .0, (21, 6) = .0, (21, 7) = .0}, datatype = float[8], order = C_order))]), ( "indepvars" ) = [x, R], ( "maxords" ) = [1, 2], ( "solvec1" ) = Vector(21, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0, (21) = .0}, datatype = float[8]), ( "startup_only" ) = false, ( "solvec4" ) = 0, ( "explicit" ) = false, ( "solmatrix" ) = Matrix(21, 7, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (1, 7) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0, (2, 7) = .0, (3, 1) = .0, (3, 2) = .0, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = .0, (3, 7) = .0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0, (4, 4) = .0, (4, 5) = .0, (4, 6) = .0, (4, 7) = .0, (5, 1) = .0, (5, 2) = .0, (5, 3) = .0, (5, 4) = .0, (5, 5) = .0, (5, 6) = .0, (5, 7) = .0, (6, 1) = .0, (6, 2) = .0, (6, 3) = .0, (6, 4) = .0, (6, 5) = .0, (6, 6) = .0, (6, 7) = .0, (7, 1) = .0, (7, 2) = .0, (7, 3) = .0, (7, 4) = .0, (7, 5) = .0, (7, 6) = .0, (7, 7) = .0, (8, 1) = .0, (8, 2) = .0, (8, 3) = .0, (8, 4) = .0, (8, 5) = .0, (8, 6) = .0, (8, 7) = .0, (9, 1) = .0, (9, 2) = .0, (9, 3) = .0, (9, 4) = .0, (9, 5) = .0, (9, 6) = .0, (9, 7) = .0, (10, 1) = .0, (10, 2) = .0, (10, 3) = .0, (10, 4) = .0, (10, 5) = .0, (10, 6) = .0, (10, 7) = .0, (11, 1) = .0, (11, 2) = .0, (11, 3) = .0, (11, 4) = .0, (11, 5) = .0, (11, 6) = .0, (11, 7) = .0, (12, 1) = .0, (12, 2) = .0, (12, 3) = .0, (12, 4) = .0, (12, 5) = .0, (12, 6) = .0, (12, 7) = .0, (13, 1) = .0, (13, 2) = .0, (13, 3) = .0, (13, 4) = .0, (13, 5) = .0, (13, 6) = .0, (13, 7) = .0, (14, 1) = .0, (14, 2) = .0, (14, 3) = .0, (14, 4) = .0, (14, 5) = .0, (14, 6) = .0, (14, 7) = .0, (15, 1) = .0, (15, 2) = .0, (15, 3) = .0, (15, 4) = .0, (15, 5) = .0, (15, 6) = .0, (15, 7) = .0, (16, 1) = .0, (16, 2) = .0, (16, 3) = .0, (16, 4) = .0, (16, 5) = .0, (16, 6) = .0, (16, 7) = .0, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = .0, (17, 5) = .0, (17, 6) = .0, (17, 7) = .0, (18, 1) = .0, (18, 2) = .0, (18, 3) = .0, (18, 4) = .0, (18, 5) = .0, (18, 6) = .0, (18, 7) = .0, (19, 1) = .0, (19, 2) = .0, (19, 3) = .0, (19, 4) = .0, (19, 5) = .0, (19, 6) = .0, (19, 7) = .0, (20, 1) = .0, (20, 2) = .0, (20, 3) = .0, (20, 4) = .0, (20, 5) = .0, (20, 6) = .0, (20, 7) = .0, (21, 1) = .0, (21, 2) = .0, (21, 3) = .0, (21, 4) = .0, (21, 5) = .0, (21, 6) = .0, (21, 7) = .0}, datatype = float[8], order = C_order), ( "depvars" ) = [theta], ( "solmat_is" ) = 0, ( "adjusted" ) = false, ( "matrixhf" ) = true, ( "norigdepvars" ) = 1, ( "stages" ) = 1, ( "theta" ) = 1/2, ( "ICS" ) = [1], ( "multidep" ) = [false, false], ( "soltimes" ) = Vector(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]), ( "depeqn" ) = [1], ( "method" ) = theta, ( "depshift" ) = [1], ( "depdords" ) = [[[1, 2]]], ( "matrixproc" ) = proc (v, vp, vpp, t, x, k, h, n, mat) local _s1, _s2, _s3, xi; _s1 := h^2; _s2 := -(3/2)/k; _s3 := (1/2)*(8*k+3*_s1)/(k*h^2); mat[3] := -(3/2)/h; mat[4] := 2/h; mat[5] := -(1/2)/h; mat[7*n-4] := 1; for xi from 2 to n-1 do mat[7*xi-4] := _s2*x[xi]^4+_s3; mat[7*xi-5] := (h-2*x[xi])/(_s1*x[xi]); mat[7*xi-3] := -(h+2*x[xi])/(_s1*x[xi]) end do end proc, ( "solution" ) = Array(1..3, 1..21, 1..1, {(1, 1, 1) = .0, (1, 2, 1) = .0, (1, 3, 1) = .0, (1, 4, 1) = .0, (1, 5, 1) = .0, (1, 6, 1) = .0, (1, 7, 1) = .0, (1, 8, 1) = .0, (1, 9, 1) = .0, (1, 10, 1) = .0, (1, 11, 1) = .0, (1, 12, 1) = .0, (1, 13, 1) = .0, (1, 14, 1) = .0, (1, 15, 1) = .0, (1, 16, 1) = .0, (1, 17, 1) = .0, (1, 18, 1) = .0, (1, 19, 1) = .0, (1, 20, 1) = .0, (1, 21, 1) = .0, (2, 1, 1) = .0, (2, 2, 1) = .0, (2, 3, 1) = .0, (2, 4, 1) = .0, (2, 5, 1) = .0, (2, 6, 1) = .0, (2, 7, 1) = .0, (2, 8, 1) = .0, (2, 9, 1) = .0, (2, 10, 1) = .0, (2, 11, 1) = .0, (2, 12, 1) = .0, (2, 13, 1) = .0, (2, 14, 1) = .0, (2, 15, 1) = .0, (2, 16, 1) = .0, (2, 17, 1) = .0, (2, 18, 1) = .0, (2, 19, 1) = .0, (2, 20, 1) = .0, (2, 21, 1) = .0, (3, 1, 1) = .0, (3, 2, 1) = .0, (3, 3, 1) = .0, (3, 4, 1) = .0, (3, 5, 1) = .0, (3, 6, 1) = .0, (3, 7, 1) = .0, (3, 8, 1) = .0, (3, 9, 1) = .0, (3, 10, 1) = .0, (3, 11, 1) = .0, (3, 12, 1) = .0, (3, 13, 1) = .0, (3, 14, 1) = .0, (3, 15, 1) = .0, (3, 16, 1) = .0, (3, 17, 1) = .0, (3, 18, 1) = .0, (3, 19, 1) = .0, (3, 20, 1) = .0, (3, 21, 1) = .0}, datatype = float[8], order = C_order), ( "totalwidth" ) = 7, ( "rightwidth" ) = 0, ( "solmat_i2" ) = 0, ( "minspcpoints" ) = 4, ( "erroraccum" ) = true, ( "eqndep" ) = [1], ( "errorest" ) = false, ( "banded" ) = true, ( "solspace" ) = Vector(21, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0, (6) = .0, (7) = .0, (8) = .0, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = .0, (14) = .0, (15) = .0, (16) = .0, (17) = .0, (18) = .0, (19) = .0, (20) = .0, (21) = 1.0}, datatype = float[8]), ( "solmat_i1" ) = 0, ( "timeadaptive" ) = false, ( "spacestep" ) = 0.500000000000000e-1, ( "initialized" ) = false, ( "vectorhf" ) = true, ( "linear" ) = true, ( "spacevar" ) = R, ( "periodic" ) = false, ( "spaceadaptive" ) = false, ( "mixed" ) = false, ( "inputargs" ) = [(-(3/2)*R^4+3/2)*(diff(theta(x, R), x))-4*(diff(theta(x, R), R)+R*(diff(diff(theta(x, R), R), R)))/R, {theta(0, R) = 1, theta(x, 1) = 0, (D[2](theta))(x, 0) = 0}], ( "bandwidth" ) = [1, 3], ( "PDEs" ) = [(-(3/2)*R^4+3/2)*(diff(theta(x, R), x))-4*(diff(theta(x, R), R)+R*(diff(diff(theta(x, R), R), R)))/R], ( "leftwidth" ) = 1 ] ); if xv = "left" then return INFO["solspace"][1] elif xv = "right" then return INFO["solspace"][INFO["spacepts"]] elif tv = "start" then return INFO["t0"] elif not (type(tv, 'numeric') and type(xv, 'numeric')) then error "non-numeric input" end if; if xv < INFO["solspace"][1] or INFO["solspace"][INFO["spacepts"]] < xv then error "requested %1 value must be in the range %2..%3", INFO["spacevar"], INFO["solspace"][1], INFO["solspace"][INFO["spacepts"]] end if; dary := Vector(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); daryt := 0; daryx := 0; dvars := []; errest := false; nd := nops(INFO["depvars"]); if dary[nd+1] <> tv then try `pdsolve/numeric/evolve_solution`(INFO, tv) catch: msg := StringTools:-FormatMessage(lastexception[2 .. -1]); if tv < INFO["t0"] then error cat("unable to compute solution for %1<%2:
", msg), INFO["timevar"], INFO["failtime"] else error cat("unable to compute solution for %1>%2:
", msg), INFO["timevar"], INFO["failtime"] end if end try end if; if dary[nd+1] <> tv or dary[nd+2] <> xv then `pdsolve/interp2dto0d`(3, INFO["soltimes"], INFO["spacepts"], INFO["solspace"], nops(INFO["depvars"]), INFO["solution"], true, tv, xv, dary); if errest then `pdsolve/interp2dto0d`(3, INFO["soltimes"], INFO["spacepts"], INFO["err_t"], nops(INFO["depvars"]), INFO["solution"], true, tv, xv, daryt); `pdsolve/interp2dto0d`(3, INFO["soltimes"], INFO["spacepts"], INFO["err_x"], nops(INFO["depvars"]), INFO["solution"], true, tv, xv, daryx) end if end if; dary[nd+1] := tv; dary[nd+2] := xv; if dvars = [] then [seq(dary[i], i = 1 .. INFO["norigdepvars"])] else vals := NULL; for i to nops(dvars) do j := eval(dvars[i]); try if errest then vals := vals, evalhf(j(tv, xv, dary, daryt, daryx)) else vals := vals, evalhf(j(tv, xv, dary)) end if catch: userinfo(5, `pdsolve/numeric`, `evalhf failure`); try if errest then vals := vals, j(tv, xv, dary, daryt, daryx) else vals := vals, j(tv, xv, dary) end if catch: vals := vals, undefined end try end try end do; [vals] end if end proc; stype := "1st"; if nargs = 1 then if args[1] = "left" then return solnproc(0, "left") elif args[1] = "right" then return solnproc(0, "right") elif args[1] = "start" then return solnproc("start", 0) else error "too few arguments to solution procedure" end if elif nargs = 2 then if stype = "1st" then tv := evalf(args[1]); xv := evalf(args[2]) else tv := evalf(args[2]); xv := evalf(args[1]) end if; if not (type(tv, 'numeric') and type(xv, 'numeric')) then if procname <> unknown then return ('procname')(args[1 .. nargs]) else ndsol := pointto(solnproc("soln_procedures")[1]); return ('ndsol')(args[1 .. nargs]) end if end if else error "incorrect arguments to solution procedure" end if; vals := solnproc(tv, xv); vals[1] end proc

(6)

NULL

gg := U(x, 1):

NULL

thm := int(U(x, R)*Uu, R = 0 .. 1):

 

 

NULL

 

Download U(R)_numériqueg2.mw

I have a nonlinear function Q(a,b,c,d,x,y) and I'd like to get the optimum (x*,y*) for different values of (a,b,c,d). The usual sintax:

NLPSolve(Q(10, 1, 5, 2, x,y), x= 0 .. 50, y = 0 .. 50, initialpoint = {x = 2,y= .5}, assume = nonnegative) does not work when Q contains numerical integration, that is evalf (Int). I have no problem with the definite integral evalf(int). The problem is that most of the cases required numerical integration so I need the former expression.

I'd appreciate very much if someone could help me.

vz := 2*(-eta^2+1);

D_im := .22;

r0 := 1;

pde := diff(vz*Y(eta, z), z)-D_im*((diff(eta*(diff(Y(eta, z), eta)), eta))/eta+diff(Y(eta, z), `$`(z, 2)))/r0 = 0;

pde := expand(%);

ibc := [Y(1, z) = 0, (D[1](Y))(0, z) = 0, Y(eta, 0) = 1, (D[2](Y))(eta, 0) = 0];

sol := pdsolve(pde, ibc, numeric, time = z, range = 0 .. 1);

pds := sol:-value(z = 0, output = listprocedure);

sol:-plot(z = 0.1e-3, numpoints = 50, color = blue, view = 0 .. 1)

So I was trying to solve this conservation equation for the radial coordinate eta and the z coordinate being treated as time. The flow is in z direction. Now unfortunately it is diverging. Not sure why though. What am I doing wrong?

Hello, My problem is as following:

 

I have tried 2 options for solving the problem below, trying to plot the behaviour of a system to a predetermined function.

First I tried to use dsolve as usual:

restart; with(plots); C := setcolors(); with(LinearAlgebra);
eq1 := Force = Mass*(diff(y(t), `$`(t, 2)));
formula1 := 2.6*BodyWeight*abs(sin(4*Pi*t));
2.6 BodyWeight |sin(4 Pi t)|
BodyWeight := 80*9.81;
plot(formula1, t = 0 .. 2);


eq2 := formula1-SpringConstant*(diff(y(t), t)) = Mass*(diff(y(t), `$`(t, 2)));
/ d \ / d /
2040.480 |sin(4 Pi t)| - SpringConstant |--- y(t)| = Mass |--- |
\ dt / \ dt \

d \\
--- y(t)||
dt //
Mass := .200;
Springt := 200;
200
SpringConstant := Youngsmodulus*Surface/DeltaLength;
DeltaLength := 0.2e-1-y(t);
Surface := .15;
Youngsmodulus := 1600*10^6-20*t^2;
eq2;
/ 2 \ / d \
0.15 \-20 t + 1600000000/ |--- y(t)|
\ dt /
2040.480 |sin(4 Pi t)| - ------------------------------------- =
0.02 - y(t)

/ d / d \\
0.200 |--- |--- y(t)||
\ dt \ dt //

incs := y(0) = 0, (D(y))(0) = 0;
eq4 := dsolve({eq2, incs});
Warning: System is inconsistent

 

Second, I tried using a numerical solving, with maxfun.


restart; with(plots); C := setcolors(); with(LinearAlgebra);
eq1 := Force = Mass*(diff(y(t), `$`(t, 2)));
formula1 := 2.6*BodyWeight*abs(sin(4*Pi*t));
2.6 BodyWeight |sin(4 Pi t)|
BodyWeight := 80*9.81;
plot(formula1, t = 0 .. 2);


eq2 := formula1-SpringConstant*(diff(y(t), t)) = Mass*(diff(y(t), `$`(t, 2)));
/ d \ / d /
2040.480 |sin(4 Pi t)| - SpringConstant |--- y(t)| = Mass |--- |
\ dt / \ dt \

d \\
--- y(t)||
dt //
Mass := .200;
Springt := 200;
200
SpringConstant := Youngsmodulus*Surface/DeltaLength;
DeltaLength := 0.2e-1-y(t);
Surface := .15;
Youngsmodulus := 1600*10^6-20*t^2;
eq2;
/ 2 \ / d \
0.15 \-20 t + 1600000000/ |--- y(t)|
\ dt /
2040.480 |sin(4 Pi t)| - ------------------------------------- =
0.02 - y(t)

/ d / d \\
0.200 |--- |--- y(t)||
\ dt \ dt //

incs := y(0) = 0, (D(y))(0) = 0;
eq4 := dsolve({eq2, incs}, y(t), type = numeric, output = listprocedure, maxfun = 10^7);
[
[t = proc(t) ... end;, y(t) = proc(t) ... end;,
[

d ]
--- y(t) = proc(t) ... end;]
dt ]

test := rhs(eq4[2]);
proc(t) ... end;

This one does plot, but no further than 0.2*10^-6. I have tried compiling the data, but this has not worked yet.

 

Does anyone know a way to work around such a problem. Is it possible to plot the equation using a for loop? If yes, how?

 

 

Hello guys,

I was just playing around with differential equations, when I noticed that symbolic solution is  different from the numerical.What is the reason for this strange behavior?


ODE := (diff(y(x), x))*(ln(y(x))+x) = 1

sol := dsolve({ODE, y(1) = 1}, y(x))

a := plot(op(2, sol), x = .75 .. 2, color = "Red");
sol2 := dsolve([ODE, y(1) = 1], numeric, range = .75 .. 2);

with(plots);
b := odeplot(sol2, .75 .. 2, thickness = 4);
display({a, b});

 

 

Strange_issue.mw

Mariusz Iwaniuk

Hello everybody.

I'm trying to obtain the numerical solution of a differential equation. Unfortunately, this prove to be quite challenging. I was able to obtain a rough solution using mathematica, but nothing more. The function is strictly increasing (for sure).

Any help is really REALLY appreciated, thanks!

 

``

deq1 := 1/(b-f(b)) = (2*(3-(1-f(b)*(diff(f(b), b, b)))/((diff(f(b), b))*(diff(f(b), b)))))/(1-2*(b-(1-f(b))/(diff(f(b), b))))

1/(b-f(b)) = 2*(3-(1-f(b)*(diff(diff(f(b), b), b)))/(diff(f(b), b))^2)/(1-2*b+2*(1-f(b))/(diff(f(b), b)))

(1)

ic1 := eval(f(b), b = 3/8) = 0, eval(f(b), b = 1/2) = 1/2

f(3/8) = 0, f(1/2) = 1/2

(2)

digits := 3

3

(3)

dsol1 := dsolve({deq1, ic1}, method = bvp[middefer], numeric, range = 3/8 .. 1/2)

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

 

``

 

Download diffeqn.mw

How can i over come convergence error, i am unable to apply approxsoln appropriately and continouation as well. regards

N := 5;

-(1/2)*Pr*n*x*(diff(f(x), x))*(diff(theta(x), x))-(1/2)*Pr*(n+1)*f(x)*(diff(theta(x), x))-(1/2)*(n+1)*(diff(diff(theta(x), x), x))+Pr*gamma*((1/4)*(n^2-3*n+3)*x^2*(diff(f(x), x))*(diff(diff(f(x), x), x))*(diff(theta(x), x))+(1/4)*(2*n^2+5*n+3)*f(x)*(diff(f(x), x))*(diff(theta(x), x))+(1/4)*n(n+1)*x*f(x)*(diff(diff(f(x), x), x))*(diff(theta(x), x))+(1/4)*(2*n^2+3*n-3)*x*(diff(f(x), x))^2*(diff(theta(x), x))+(1/4)*(n-1)*x^2*(diff(diff(f(x), x), x))*(diff(theta(x), x))+(1/2)*n*(n+1)*x*f(x)*(diff(f(x), x))*(diff(diff(theta(x), x), x))+(1/4)*(n^2-1)*(diff(f(x), x))^2*(diff(theta(x), x))+(1/4)*(n+1)^2*f(x)^2*(diff(diff(theta(x), x), x))+(1/4)*(n-1)^2*x^2*(diff(f(x), x))^2*(diff(diff(theta(x), x), x))) = 0

(1)

bc := (D(theta))(0) = -Bi*(1-theta(0)), theta(N) = 0, f(0) = 0, (D(f))(0) = 0, (D(f))(N) = 1;

(D(theta))(0) = -Bi*(1-theta(0)), theta(5) = 0, f(0) = 0, (D(f))(0) = 0, (D(f))(5) = 1

(2)

a1 := dsolve(subs(beta = .1, n = .5, Pr = 10, gamma = .1, Bi = 50, {bc, eq1, eq2}), numeric, method = bvp[midrich], abserr = 10^(-8), output = array([seq(.1*i, i = 0 .. 10*N)]))

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

 

``

 

Download ehtasham.mwehtasham.mw

i have attcahed my ode with complex bvp

can anyone solved mine

NULL

restart

with(plots):

NULL

Eq1 := (11-10*d)*(diff(h(eta), eta))+2*f(eta) = 0;

(11-10*d)*(diff(h(eta), eta))+2*f(eta) = 0

 

(11-10*d)*(diff(diff(f(eta), eta), eta))-h(eta)*(diff(f(eta), eta))-f(eta)^2+g(eta)^2 = 0

 

diff(diff(g(eta), eta), eta)-h(eta)*(diff(g(eta), eta))-2*f(eta)*g(eta) = 0

 

diff(p(eta), eta)+2*(diff(f(eta), eta))-2*f(eta)*h(eta) = 0

(1)

NULL

NULL

`V&lambda;` := [0.5e-1, 1.5, 1.5]:

etainf := 3:

bcs := h(0) = 0, p(0) = 0, (D(f))(0) = lambda*f(0)^(4/3)/(f(0)^2+(1-g(0))^2)^(1/3), (D(g))(0) = -Typesetting:-delayDotProduct(lambda*f(0)^(1/3)*(1-g(0)), 1/(f(0)^2+(1-g(0))^2)^(1/3)), f(etainf) = 0, g(etainf) = 0;

h(0) = 0, p(0) = 0, (D(f))(0) = lambda*f(0)^(4/3)/(f(0)^2+(1-g(0))^2)^(1/3), (D(g))(0) = -f(0)^(1/3)*(1-g(0))*lambda/(f(0)^2+(1-g(0))^2)^(1/3), f(3) = 0, g(3) = 0

(2)

NULL

dsys := {Eq1, Eq2, Eq3, Eq4, bcs}:

for i to 3 do lambda := `V&lambda;`[i]; dsol[i] := dsolve(dsys, numeric, continuation = d); print(lambda); print(dsol[i](0)) end do

Error, (in dsolve/numeric/bvp) singularity encountered

 

NULL

NULL

NULL

 

Download compre1.mw

and attch back

1 2 3 4 5 6 7 Last Page 1 of 22