Items tagged with ode

Feed

hello everyone. i wanna ask. how to do a coding of stability region on diagonally implicit Runge-Kutta for order 3 and order 4? Thank you in advanced :) 

i want to solve this equation,

y''(x)=5*exp(-10/y'(x)) on ]0,15[ with y(0)=0,y(15)=2 

can any one help me ? thank you

Hi, I'm trying to solve this ode:
restart; with(plots); with(DEtools);

l := t -> 0.5*tanh(0.5*t);

deq := diff(f(t), t)*l(t)*(diff(f(t), t, t)*l(t)+9.8*sin(f(t)))+diff(l(t), t)*(diff(f(t), t)^2*l(t)-9.8*cos(f(t))+4*(l(t)-0.5)) = 0;

sol := dsolve({deq, f(0) = 0, D(f)(0) = 0.1}, f(t), numeric);

 

but getting an error:

Error, (in dsolve/numeric/checksing) ode system has a removable singularity at t=0. Initial data is restricted to {f(t) = 1.77632183122019}
 

How can I possibly fix this?

Respected member!

Please help me in finding the solution of this problem....
 

NULL

 

 

NULL

>   

``

NULL

restart

with(RealDomain):

r := .2:

k := 5;

5

(1)

BCSforNum1 := u(0) = 0, (D(u))(0) = 1+beta*(((D@@2)(u))(0)-(D(u))(0)*RealDomain:-`^`(k, -1)), (D(u))(m) = 0, ((D@@2)(u))(m) = 0;

u(0) = 0, (D(u))(0) = 1+.2*((D@@2)(u))(0)-0.4000000000e-1*(D(u))(0), (D(u))(6) = 0, ((D@@2)(u))(6) = 0

 

v(0) = 1, v(6) = 0

(2)

numsol1 := dsolve({BCSforNum1, BCSforNum2, ODEforNum1, ODEforNum2}, numeric, output = listprocedure)

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

 

``


 

Download mplprimes.mw

Dear please check once it showing an error program.mw as intial value is not conververging

Dear sir,

I tried to solve a fourth order problem. But I got the error message as better to use midpoint method. Can I know what is midpoint method and here I uploading the problem please verify it if I did anything mistake?program.mw

i am interested to numerically solve the 3 non-linear coupled ODE's in the 3 different intervals of rho(define in attached file) for the different corresponding parameters alpha and beta at i=0..2.

density_interval_solution.mw 

Download 2222222222222222222.mw

 

> restart; 
> A[0] := 10^(-3); a := 10^5;
> sys := diff(R(theta), theta) = A[0]*exp(2*mu(theta))*sin(theta)/(2*a), R(theta) = 2*exp(-2*mu(theta))*(1-(diff(mu(theta), `$`(theta, 2)))-cot(theta)*(diff(mu(theta), theta)));
> cond := R(0) = 10^(-5), mu(0) = 118.92, (D(mu))(0) = 0;
> F := dsolve({cond, sys}, [R(theta), mu(theta)], numeric);
> with(plots);
> odeplot(F, [theta, R(theta)], 0 .. 3.14, color = black, thickness = 3, linestyle = 4)
> odeplot(F, [theta, mu(theta)], 0 .. 3.14, color = blue, thickness = 3, linestyle = 1)

After last two lines maple writes:

Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

And gives me empty plots. I can't figure out where an error can be. Some things I noticed:

Maple doesn't calculate the system before and after zero. If I change the range from 0..3.14 to -10..10 or to 0.00001..0.00001, it gives me 2 errors for 1 plot.

Also if I change the condition mu(0) = 118.92 to mu(0) = 1 or mu(0) = 50 or mu(0) = 80, it works. After ~80 it gives an error. I can't imagine where could appear a division by 0 or some other mistake.

Dear All,

The following code calculates the Rotation_number for given set of parameters. I would like to plot the value of Rotation_number against R for R from 500 to 1000. Since there are many variables for each iteration, I would like to know if there is a way to create a for loop for the following code while automatically clear variables except for Rotation_number to speed it up? 


 

``

restart:
with(plots): with(DEtools): with(plottools):with(LinearAlgebra):
t_start:=2500: t_end:=5000:

b:=-3: c:=3: v1:=1: f:=-4: v2:=2.0: omega:=1: epsilon:=evalf(1/R): R:=100: k:=0.25:

kH:=(c+1)/(b-1):

sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u1(t),
     diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2-b*z(t)^2)*u2(t),
     diff(z(t),t)=z(t)*(kH*v1+c*u1(t)^2+c*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):

solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},
             {u1(t),u2(t),z(t)},
              numeric, method=rkf45, maxfun=10^7,
              events=[[[u1(t)=0, u2(t)>0], halt]]):
evs:=Array():
evs(1,1..4):=Array([t,u1(t),u2(t),z(t)]):
interface(warnlevel=0):

for k from 2 do
    w:=solA(t_end):
    if rhs(w[1])<t_end
    then evs(k,1..4):=Array(map(rhs, w));
         solA(eventclear);
    else break;
    fi
od:

interface(warnlevel=3):
M:=DeleteRow(convert(evs,matrix),1):

T:=M[..,1]:

for i from 1 do
tt:=T[i]:
if tt>=t_start then break; end if;
end do:

R:=M[i..,3]: Z:=M[i..,4]:

N_alpha:=numelems(R):

r_center:=add(R[i],i=1..N_alpha)/N_alpha:
z_center:=add(Z[i],i=1..N_alpha)/N_alpha:

Zshift:=Z-~r_center: Rshift:=R-~z_center:

alpha:=(arctan~(Zshift/~Rshift)+(Pi/2)*sign~(Rshift))/~(2*Pi)+~0.5:



alpha2:=alpha[2..N_alpha]:
alpha1:=alpha[1..N_alpha-1]:

del_alpha:=alpha2-alpha1:
m:=numelems(del_alpha):

for i from 1 to m do

if del_alpha[i]<0
then del_alpha[i]:=del_alpha[i]+1;
else del_alpha[i]:=del_alpha[i];
fi:
od:

Rotation:=add(del_alpha[i],i=1..m)/m;

HFloat(0.8000000001553702)

(1)

NULL

Rotation_number.mw

Thank you!

Very kind wishes,

Wang Zhe

TQ.mw

Can any one help for finding the solution of these differntial equations and then plotting the graph for differnt values of M

(FILE ATTACHED)
 

eqn1 := (R/(R-theta(eta))+Omega)*(diff(f(eta), `$`(eta, 3)))+f(eta)*(diff(f(eta), `$`(eta, 2)))+R*(diff(f(eta), `$`(eta, 2)))*(diff(theta(eta), eta))/(R-theta(eta))^2+Omega*(diff(g(eta), eta))+lambda*theta(eta)*Cos(alpha)-M*(diff(f(eta), eta))^2 = 0; eqn2 := (R/(R-theta(eta))+(1/2)*Omega)*(diff(g(eta), `$`(eta, 2)))-2*Omega*(2*g(eta)+diff(f(eta), `$`(eta, 2)))+(diff(f(eta), eta))*g(eta)+(diff(g(eta), eta))*f(eta)+R*(diff(g(eta), eta))*(diff(theta(eta), eta))/(R-theta(eta))^2 = 0; eqn3 := (1+`&epsilon;`*theta(eta))*(diff(theta(eta), `$`(eta, 2)))+`&epsilon;`*(diff(theta(eta), eta))^2+Pr*(f(eta)*(diff(theta(eta), eta))-(diff(f(eta), eta))*theta(eta))+Q*theta(eta)+L*exp(-eta) = 0

(R/(R-theta(eta))+Omega)*(diff(diff(diff(f(eta), eta), eta), eta))+f(eta)*(diff(diff(f(eta), eta), eta))+R*(diff(diff(f(eta), eta), eta))*(diff(theta(eta), eta))/(R-theta(eta))^2+Omega*(diff(g(eta), eta))+lambda*theta(eta)*Cos(alpha)-M*(diff(f(eta), eta))^2 = 0

 

(R/(R-theta(eta))+(1/2)*Omega)*(diff(diff(g(eta), eta), eta))-2*Omega*(2*g(eta)+diff(diff(f(eta), eta), eta))+(diff(f(eta), eta))*g(eta)+(diff(g(eta), eta))*f(eta)+R*(diff(g(eta), eta))*(diff(theta(eta), eta))/(R-theta(eta))^2 = 0

 

(1+epsilon*theta(eta))*(diff(diff(theta(eta), eta), eta))+epsilon*(diff(theta(eta), eta))^2+Pr*(f(eta)*(diff(theta(eta), eta))-(diff(f(eta), eta))*theta(eta))+Q*theta(eta)+L*exp(-eta) = 0

(1)

Omega := 2.; M := .5; R := 5; lambda := 20; `&epsilon;` := .2; Pr := 1; Q := .5; L := .5; W := .5; n := .1; alpha := (1/6)*Pi

bc := f(0) = W, (D(f))(0) = 0, (D(f))(infinity) = 0, (D(theta))(0) = -1, theta(infinity) = 0, g(0) = -n*(DD(f))(0), g(infinity) = 0

f(0) = W, (D(f))(0) = 0, (D(f))(N) = 0, (D(theta))(0) = -1, theta(N) = 0, g(0) = -n*(DD(f))(0), g(N) = 0

(2)

``


 

Download TQ.mw

 

hello agian i have the following problem. i need to fit data using the following model:

ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;

ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;

ode_P2 := diff(P2(t), t) = -k2*(-keq*P1(t)+P2(t))/keq-k4*P2(t)-P2(t)/T1_p2;

ode_P2e := diff(P2_e(t), t) = k4*P2(t)-P2_e(t)/T1_p2_e;

ode_system := ode_sub, ode_P1, ode_P2, ode_P2e

known paramters: s0 := 10000; k2 := 1000; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1

initial conditions: init := S(0) = s0, P1(0) = 0, P2(0) = 0, P2_e(0) = 0

using the following data the fitting is fine:

T := [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100]

data_s := [10000.00001377746, 7880.718836217512, 6210.572917625112, 4894.377814704496, 3857.121616806618, 3039.689036293312, 2395.493468824973, 1887.821030424934, 1487.738721378254, 1172.445015238064, 923.970970533911, 728.1555148262845, 573.8389058920359, 452.2263410725434, 356.3868133709972, 280.8584641247961, 221.3366863178145, 174.4291648589732, 137.4627967029932, 108.3305246342751, 85.37230370803576, 67.2794974831867, 53.0210755540487, 41.78436556408739, 32.92915589156605, 25.95049906181449, 20.45092590987601, 16.11678944747619, 12.70115104646253, 10.009439492356, 7.888058666357939, 6.216464754698855, 4.899036441885205, 3.860793948052506, 3.042584031379331, 2.397778731385364, 1.889601342133927, 1.489145259940784, 1.173591417634974, .9248694992255094, .7288443588090404, .5743879613705721, .4526024635132348, .3567687570029392, .2810508404011898, .2215650452813882, .174603181767829, .1376243372528356, .1084232889842753, 0.8542884707952822e-1, 0.6738282660463157e-1]

data_p1 := [0.1194335334401124e-4, 244.8746046039949, 374.8721199398692, 430.5392805383767, 439.6598364813143, 421.0353424914179, 387.1842556288343, 346.2646897222593, 303.4377508746471, 261.8283447091155, 223.1996547160051, 188.4213144493491, 157.7924449350029, 131.2622073983344, 108.5771635112278, 89.37951009190863, 73.26979150957087, 59.84578653950572, 48.72563358658898, 39.56010490461378, 32.03855466968536, 25.88922670933322, 20.87834763772145, 16.80708274458702, 13.50774122768974, 10.84014148654258, 8.687656394505874, 6.954093898245485, 5.560224107929433, 4.441209458524726, 3.544128529596104, 2.825755811619965, 2.251247757181308, 1.792233305651086, 1.425861347838012, 1.133566009019768, .90081361320016, .7153496336919163, .5678861241754847, .4505952916932289, .3572989037753538, .2832489239941939, .2244289868248577, .1778450590752305, .1408633578784151, .1114667192753896, 0.8826814044702111e-1, 0.6979315954603076e-1, 0.5526502783606788e-1, 0.4370354298880999e-1, 0.3456334307662573e-1]

data_p2_p2e := [-0.1821397630630296e-4, 1000.40572909871, 1568.064904416198, 1848.900129881268, 1944.147939710225, 1923.352973299286, 1833.705314342611, 1706.726235937363, 1563.036042115902, 1415.741121363331, 1272.825952816517, 1138.833091575137, 1016.03557293539, 905.2470623752856, 806.3754051707843, 718.7979384094563, 641.6091822001032, 573.7825966275556, 514.2718125966452, 462.0710416682647, 416.2491499591012, 375.9656328260581, 340.4748518743264, 309.124348579787, 281.3484612079911, 256.6599441255977, 234.6416876403397, 214.9378461256197, 197.2453333305823, 181.3059798685686, 166.90059218783, 153.8422174795751, 141.9717783871606, 131.1529759058513, 121.2692959115991, 112.2199678247825, 103.9187616370328, 96.29007054510998, 89.26877137303566, 82.79743967408479, 76.82586273439793, 71.30944943617081, 66.2085909515396, 61.48838150744805, 57.11714763242225, 53.06666224006544, 49.31130738219119, 45.82807853990728, 42.59597194910467, 39.59575450632008, 36.81013335261527]

the fitting is done the following way: 

P1fu, P2fu, P2e_fu, Sfu := op(subs(res, [P1(t), P2(t), P2_e(t), S(t)]))

making residuals:

Q := proc (T1_p1, k1, keq, k4) local P1v, P2v, P2e_v, Sv, resid; option remember;

res(parameters = [T1_p1, k1, keq, k4]);

try P1v := `~`[P1fu](T); P2v := `~`[P2fu](T); P2e_v := `~`[P2e_fu](T); Sv := `~`[Sfu](T); resid := [P1v-data_p1, P2v+P2e_v-data_p2_p2e, Sv-data_s]; return [seq(seq(resid[i][j], i = 1 .. 3), j = 1 .. nops(T))] catch: return [1000000$3*nops(T)] end try end proc;
q := [seq(subs(_nn = n, proc (T1_p1, k1, keq, k4) Q(args)[_nn] end proc), n = 1 .. 3*nops(T))];

finding inital point for the LSsolve:

L := fsolve(q[2 .. 5], [10, 0.2e-1, 4, 4])

fitting the data with the intial point: 

solLS := Optimization:-LSSolve(q, initialpoint = L)

this is all good however, when i used the following data it did not turn out so well (using the same approch as above):

data_s := [96304.74567, 77385.03700, 62621.83067, 51239.94333, 42663.82367, 35084.74100, 28480.28367, 23066.01467, 18774.73700, 15179.13700, 12278.50767, 9937.652000, 8046.848333, 6521.242000, 5287.811667, 4277.779000, 3466.518333, 2835.467000, 2297.796333, 1861.249667, 1529.654000, 1235.353000, 999.6626667, 826.2343333, 667.9480000, 559.9230000, 449.2790000, 376.4860000, 289.1203333, 236.1483333]

data_p1 := [0.86e-1, 3.904, 26.975, 31.719, 41.067, 46.779, 52.115, 43.101, 44.344, 41.094, 36.523, 27.742, 26.543, 28.062, 22.178, 21.303, 14.951, 17.871, 11.422, 12.051, 9.232, 6.817, 6.1, .717, 1.215, 6.146, .772, .375, 2.595, .518]

data_p2_p2e := [-3.024, 22.238, 61.731, 103.816, 132.695, 159.069, 167.302, 160.188, 158.398, 152.943, 146.745, 135.22, 132.145, 120.413, 107.864, 95.339, 90.775, 81.828, 71.065, 70.475, 62.872, 49.955, 40.858, 42.938, 41.311, 35.583, 31.573, 29.841, 29.558, 21.762]

the known parameters is the case are:

s0 := 96304.74567; k2 := 10^5; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1

additionally the following fuction affects the solution: 

K:=t->cos((1/180)*beta*Pi)^(t/Tr)

i included this by doing this:

P1fu_K := proc (t) options operator, arrow; P1fu(t)*K(t) end proc;

P2fu_K := proc (t) options operator, arrow; P2fu(t)*K(t) end proc;

P2e_fu_K := proc (t) options operator, arrow;

P2e_fu(t)*K(t) end proc;

Sfu_K := proc (t) options operator, arrow; Sfu(t)*K(t) end proc

resulting in the following residuals:

Q := proc (T1_p1, k1, keq, k4) local P1v, P2v, P2e_v, Sv, resid; option remember;

res(parameters = [T1_p1, k1, keq, k4]);

try P1v := `~`[P1fu_K](T); P2v := `~`[P2fu_K](T); P2e_v := `~`[P2e_fu_K](T); Sv := `~`[Sfu_K](T); resid := [P1v-data_p1, P2v+P2e_v-data_p2_p2e, Sv-data_s]; return [seq(seq(resid[i][j], i = 1 .. 3), j = 1 .. nops(T))] catch: return [1000000$3*nops(T)] end try end proc;
q := [seq(subs(_nn = n, proc (T1_p1, k1, keq, k4) Q(args)[_nn] end proc), n = 1 .. 3*nops(T))];

i think my problem is that the inital point in this case is not known. all i know is that all the fitted parameters should be positive and that k1<1, k4>10, keq>1 and T1_p1>100 (more i do not know) - is there a way to determin the inital point without guessing?

i also know the results, which should be close to these values:

k1=0.000438, k4=0.0385, keq=2.7385 and T1_p1=36.8 the output fit should look something like this

where the red curve is (PP2(t)+PP2_e(t))*K(t)

the blue cure is PP1(t)*K(t) 

anyone able to help - i've tried for 2 days now. it might be that  ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1 should be changed into ode_P1 := diff(P1(t), t) = k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;

i've tried this but it didnt seem to do much 

anyone able to help?:)

NB stiff=true can be used within the dsolve to speed up the process if needed:)

 

 

 

 


 

 

restart;

Digits := 18;
with(LinearAlgebra);
f := proc (n) 3*sin(x[n]) end proc;

g := proc (n) 3*cos(x[n])

end proc;

#problem call.
for n from 0 to 0 do

e1 := expand(-y[n+3/2]+y[n]-3*y[n+1/2]+3*y[n+1]+1/11612160*(5856*h^4*g(n+1/2)-19968*h^4*g(n+3/2)+2343*h^4*g(n)-76356*h^4*g(n+1)-7058*h^4*g(n+2)+608864*h^3*f(n+1/2)+104864*h^3*f(n+3/2)+28489*h^3*f(n)+702864*h^3*f(n+1)+6439*h^3*f(n+2)));

e2 := expand(-y[n+2]+3*y[n]-8*y[n+1/2]+6*y[n+1]+1/5806080*(18768*h^4*g(n+1/2)-32880*h^4*g(n+3/2)+3867*h^4*g(n)-76356*h^4*g(n+1)-2229*h^4*g(n+2)+965728*h^3*f(n+1/2)+461728*h^3*f(n+3/2)+45953*h^3*f(n)+1405728*h^3*f(n+1)+23903*h^3*f(n+2)));

e3 := expand(-z[n]+(1/383201280*(-4207440*h^4*g(n+1/2)-930192*h^4*g(n+3/2)+371973*h^4*g(n)-3631932*h^4*g(n+1)-41259*h^4*g(n+2)+16136096*h^3*f(n+1/2)+3866720*h^3*f(n+3/2)+5752543*h^3*f(n)+5810400*h^3*f(n+1)+367681*h^3*f(n+2))+4*y[n+1/2]-3*y[n]+y[n+1])/h);

e4 := expand(-z[n+1/2]+(1/191600640*(376320*h^4*g(n+1/2)+118896*h^4*g(n+3/2)-29469*h^4*g(n)+532764*h^4*g(n+1)+5079*h^4*g(n+2)-5812112*h^3*f(n+1/2)-508016*h^3*f(n+3/2)-381553*h^3*f(n)-1236168*h^3*f(n+1)-45511*h^3*f(n+2))-y[n]+y[n+1])/h);

e5 := expand(-z[n+1]+(1/383201280*(-31920*h^4*g(n+1/2)-433776*h^4*g(n+3/2)+71547*h^4*g(n)-2519748*h^4*g(n+1)-17493*h^4*g(n+2)+18565216*h^3*f(n+1/2)+1933216*h^3*f(n+3/2)+885665*h^3*f(n)+10391328*h^3*f(n+1)+158015*h^3*f(n+2))-5*y[n+1/2]+y[n]+3*y[n+1])/h);

e6 := expand(-z[n+3/2]+(1/95800320*(250224*h^4*g(n+1/2)-730680*h^4*g(n+3/2)+61266*h^4*g(n)-1526256*h^4*g(n+1)-22044*h^4*g(n+2)+15680504*h^3*f(n+1/2)+4712456*h^3*f(n+3/2)+735469*h^3*f(n)+22576428*h^3*f(n+1)+203623*h^3*f(n+2))-8*y[n+1/2]+3*y[n]+5*y[n+1])/h);

e7 := expand(-z[n+2]+(1/383201280*(3873264*h^4*g(n+1/2)+332976*h^4*g(n+3/2)+497649*h^4*g(n)-1407564*h^4*g(n+1)-720255*h^4*g(n+2)+114710816*h^3*f(n+1/2)+93716192*h^3*f(n+3/2)+5705827*h^3*f(n)+191366496*h^3*f(n+1)+9635389*h^3*f(n+2))-12*y[n+1/2]+5*y[n]+7*y[n+1])/h);

e8 := expand(-p[n]+(1/191600640*(13423440*h^4*g(n+1/2)+3068304*h^4*g(n+3/2)-1621317*h^4*g(n)+11615292*h^4*g(n+1)+137451*h^4*g(n+2)-32503712*h^3*f(n+1/2)-12664928*h^3*f(n+3/2)-32539039*h^3*f(n)-16869600*h^3*f(n+1)-1223041*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e9 := expand(-p[n+1/2]+(1/191600640*(-3053856*h^4*g(n+1/2)-213216*h^4*g(n+3/2)+98049*h^4*g(n)-509436*h^4*g(n+1)-10191*h^4*g(n+2)-1045120*h^3*f(n+1/2)+831104*h^3*f(n+3/2)+1331083*h^3*f(n)-1207008*h^3*f(n+1)+89941*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e10 := expand(-p[n+1]+(1/63866880*(194160*h^4*g(n+1/2)-373968*h^4*g(n+3/2)+52329*h^4*g(n)-2514924*h^4*g(n+1)-14727*h^4*g(n+2)+14006304*h^3*f(n+1/2)+1695712*h^3*f(n+3/2)+634955*h^3*f(n)+15463008*h^3*f(n+1)+133461*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e11 := expand(-p[n+3/2]+(1/191600640*(1491168*h^4*g(n+1/2)-4758240*h^4*g(n+3/2)+190977*h^4*g(n)-509436*h^4*g(n+1)-103119*h^4*g(n+2)+46274944*h^3*f(n+1/2)+48151168*h^3*f(n+3/2)+2215307*h^3*f(n)+93985056*h^3*f(n+1)+974165*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2);

e12 := expand(-p[n+2]+(1/191600640*(4772688*h^4*g(n+1/2)+11719056*h^4*g(n+3/2)+338619*h^4*g(n)+11615292*h^4*g(n+1)-1822485*h^4*g(n+2)+59770976*h^3*f(n+1/2)+79609760*h^3*f(n+3/2)+3528289*h^3*f(n)+109647648*h^3*f(n+1)+34844287*h^3*f(n+2))-8*y[n+1/2]+4*y[n]+4*y[n+1])/h^2) end do;
M := {e || (1 .. 12)};

y_init := 1;

z_init := 0;

p_init := -2;

x_init := 0; A := 0; B := 1; N := 40;

h := evalf((B-A)/N); count := 1;

X := y[k], y[k+1/2], y[k+1], y[k+3/2], z[k], z[k+1/2], z[k+1], z[k+3/2], p[k], p[k+1/2], p[k+1], p[k+3/2];

step := seq(eval(x, x = n*h), n = 1 .. N);

y_exact := ([seq])(eval(3*cos(x)+(1/2)*x^2-2, x = n*h), n = 1 .. N);

z_exact := ([seq])(eval((1/3*(3*x^2+6*x+3))/(x^3+3*x^2+3*x+1), x = n*h), n = 1 .. N);

p_exact := ([seq])(eval((1/3*(6*x+6))/(x^3+3*x^2+3*x+1)-(1/3)*(3*x^2+6*x+3)^2/(x^3+3*x^2+3*x+1)^2, x = n*h), n = 1 .. N);
vars := seq(X, k = 1);
printf("\n%4s%13s%15s%15s\n", "@", "y_Num", "y_Exact", "y_Error");

for q to N do

for ix to 4 do

x[ix] := h*ix+x_init end do;

result := eval(`<,>`(vars), fsolve(eval(M, [x[0] = x_init, x[1/2] = x_init, x[3/2] = x_init, y[0] = y_init, y[1/2] = y_init, y[3/2] = y_init, z[0] = z_init, z[1/2] = z_init, z[3/2] = z_init, p[0] = p_init, p[1/2] = p_init, p[3/2] = p_init]), {vars}));

for k to 4 do

printf("%5.2f %14.15f", step[count], result[k]);

printf("%20.15f %10.18G \n", y_exact[count], abs(result[k]-y_exact[count]));

count := count+1;

P := [result[k]]

end do;

x_init := x[ix-1];

y_init := result[4];

z_init := result[8];

p_init := result[12]

end do;

 

please that is the code i write to solve the problem after using the matrix form to generate the value but is given me error of the form


   @        y_Num        y_Exact        y_Error
Error, invalid input: eval received fsolve({-6398.00004614630940+6400.00000000000000*y[1], -6397.99992849910140+6400.00000000000000*y[1], -6397.99909739580050+6400.00000000000000*y[1], -199.999989717789185+200.000000000000000*y[1], -40.0000000791700798+40.0000000000000000*y[1], -2.99999993737911015+3*y[1], 39.999999768462113+40.0000000000000000*y[1], -p[1]-6399.99961623646730+6400.00000000000000*y[1], -p[2]-6399.99798489466010+6400.00000000000000*y[1], -y[2]-4.99999972458202552+6*y[1], -z[1]-159.999999048856193+120.000000000000000*y[1], -z[2]-279.999973921987948+280.000000000000000*y[1]}, {p[1], p[2], p[3/2], p[5/2], y[1], y[2], y[3/2], y[5/...
 

Hi,

I need your help to classify the follwing set {0}, {1} and [0,1] are local attractor or not and in the case of local attractor how can we determine the bassin of attraction. 

ode:=diff(x(t),t)=sqrt(x(t));

how can we prove using maple which of {0}, {1} and [0,1] are local attarctor or not.

Many thanks

 

Dear All,

I would like to plot the probability density function of a state variable obtained from solving differential equations. I have found that there are functions called "PDF" and "KernelDensityPlot" in the Statistics package, but they really confuse me. Could you please point me out? My code is as follows.

Ps. Is it possible to plot the PDF directly from the solution of dsolve() without discretizing the results? 


restart:
with(plots): with(DEtools): with(plottools):with(LinearAlgebra): with(Statistics):

v1:=1: f:=-4: v2:=2.515: omega:=1: epsilon:=0.001: k:=0:
sys:=diff(u1(t),t)=v1*u1(t)-(omega+k*u2(t)^2)*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u1(t),
     diff(u2(t),t)=(omega+k*u1(t)^2)*u1(t)+v1*u2(t)-(u1(t)^2+u2(t)^2+3*z(t)^2)*u2(t),
     diff(z(t),t)=z(t)*(-v1+3*u1(t)^2+3*u2(t)^2+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4):

t_start:=50: t_end:=300: dt:=0.05: fs:=1/dt:

solA:=dsolve({sys, u1(0)=0.6, u2(0)=0.6, z(0)=0.1},
             {u1(t),u2(t),z(t)},
              type=numeric, method=rkf45, maxfun=0,
              output=Array([seq(i,i=t_start..t_end, dt)])):

u1:=solA[2,1][..,2]:
u2:=solA[2,1][..,3]:
z:=solA[2,1][..,4]:

u0:=sqrt~(u1^~2+u2^~2+z^~2):

Phi:=z/~u0:

# How could I plot the probability density of Phi (y-axis) against Phi(x-axis)?

Probability_density_function_plot.mw

Thank you!

Very kind wishes,

Wang Zhe

Respected member!
Please help me to find the solution of attached problem.

 


> subject to boundary conditions


``

 

 

NULL

restart

alpha := evalf(2*Pi*(1/180)); EP := .2; lambda := .1; HA := 5; RE := 20

ODEforNum := (1+EP)*(((D@@3)(F))(r)+4*alpha^2*(D(F))(r))+2*alpha*RE*F(r)*(D(F))(r)-HA*alpha^2*(D(F))(r)-3*EP*lambda*((1/2)*(D(F))(r)^2*((D@@3)(F))(r)+(D(F))(r)*((D@@2)(F))(r)^2)/alpha^2-EP*lambda*(72*F(r)^2*(D(F))(r)+2*(D(F))(r)^3+32*F(r)*(D(F))(r)*((D@@2)(F))(r)+2*F(r)^2*((D@@3)(F))(r)) = 0

1.2*((D@@3)(F))(r)-0.243693936e-3*(D(F))(r)+1.396263402*F(r)*(D(F))(r)-24.62104762*(D(F))(r)^2*((D@@3)(F))(r)-49.24209525*(D(F))(r)*((D@@2)(F))(r)^2-1.44*F(r)^2*(D(F))(r)-0.4e-1*(D(F))(r)^3-.64*F(r)*(D(F))(r)*((D@@2)(F))(r)-0.4e-1*F(r)^2*((D@@3)(F))(r) = 0

 
 

 

BCSforNum := F(0) = 1, (D(F))(0) = 0, F(1) = 0

Digits := 15

15

(2)

numsol := dsolve({BCSforNum, ODEforNum}, numeric, output = listprocedure)

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

 
 

 

 

 

 

 

 

 

Download MapleN.mw

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