Items tagged with ode ode Tagged Items Feed

This is the code  hw2_final.mw

 

Let me explain it.

I am sure that the mistakes must be about the expresstion of the I1(t) and I2(t). Actually if you delete I1(t) and I2(t) , the whole code works and get the picture at the bottom. 

What I want is to put the expresstion of the I1(t) and I2(t) into 'sol:=dsolve...' and 'plots...' to get the picture of I1(t) and I2(t) with respect to t. Before the t* which subject to Phi(t*)=0 (The blue line in the picture at the bottom is Phi) I want I1(t) and after t* I want I2(t).

I1(t) = (int(sqrt(2*(H(t)+omega*cos(q(t)))), q(t) = q(t)-2*Pi .. q(t), numeric))/Pi.    what I want of this experesstion is to get  'int(sqrt(2*(H(t)+omega*cos(q(t)))' from  'q(t)-2*Pi' to 'q(t)' by numeric method.This q(t) is the solution of the ODE sys.

For example(the number I used is not true,just for example) , at the point t=20, q(t)=30-2*Pi.

so I1(t)= (int(sqrt(2*(H(t)+omega*cos(x))), x = 30-2*Pi .. 30, numeric))/Pi.The I2(t) I want is similar to I1(t).

 

How can I solve it?

I'm trying to solve the differential equation.

Eq := diff(y(x), x, x) = -(x^2+1)*y(x)+K;

dsolve({Eq, y(-1) = 0, y(1) = 0}, y(x));

But this not work very well.

Best Regards,

I want to solve an ODE from Game Theory, the Cournot competition.

It says

p(q1+r2(q1))+p'(q1+r2(q1))*r2(q1)-c2'(r2(q1))=0

 where, I think,

' means diff(,q1),

c2(q2)=c*q2 for a fixed c in [0,1]

and

p(q)=max(0,1-q).

So c2,p and r2 are functions.r2 goes from [0,inf) to [0,inf).

I look for r2, which should be r2(q1)=(1-q1-c)/2 when correctly solved.

However, the command dsolve says Error in dsolve (divison by 0).

 What is wrong? How do I obtain the solution for r2 in Maple?

 

My goal is to plot the integral J with respect to t and as you can see J is a piecewise function.

This is my code.   hw2_numerical_2.mw

 

Actually it's a problem about adiabatic invariants.

If you want to know the backgroud please see this link.

www.mapleprimes.com/questions/206645-How-To--Numerically--Solve-It-

Hello,

I'm writing a code and I seem to have an issue when trying to implement a procedure. Here is the code:

 

with(plots):

Z := 75; A := 189; k := 14.6; Rm := 8*R; r0 := 10^(-8)*R; c := 137.036; ms := 105.66/(.51100)

fmtoau := 10^(-15)/(0.529177e-1*10^(-9)):

R := 1.1*fmtoau*A^(1.0/(3.0)):

f := proc (x) options operator, arrow; 1/(1+exp(k*(x-1))) end proc:

n0 := 3*Z*k^2/(4*Pi*(Pi^2+k^2)*R^3):

n := proc (r) options operator, arrow; 4*Pi*n0*f(r/R) end proc:

int(r^2*n(r), r = 0 .. Rm);

73.00000007

(1)

plot(n/n0, 0 .. 2*R);

 

v1 := unapply(int(x^2*f(x), x), x):

Vfermi := proc (r) options operator, arrow; -4*Pi*n0*R^2*(R*(v1(r/R)-v10)/r+v2Rm-v2(r/R)) end proc:

Vuniform := proc (r) options operator, arrow; piecewise(r < R, -Z*(3/2-(1/2)*r^2/R^2)/R, -Z/r) end proc:

plot([Vuniform(r), Vfermi(r), -Z/r], r = r0 .. 2*R, V = 1.2*Vfermi(r0) .. 0, legend = ["uniform charge", "Fermi distribution", "point charge"]);

 

``

plotsol1s := proc (E, K, r0, S, col) local Eqns, ICs, fnl, gnl, r, soln; global ms, c; Eqns := diff(fnl(r), r) = gnl(r)*[E+2*ms*c^2-Vuniform(r)]/c-(K+1)*fnl(r)/r, diff(gnl(r), r) = -fnl(r)*[E-Vuniform(r)]/c-(1-K)*fnl(r)/r; ICs := fnl(r0) = 1, gnl(r0) = 0; soln := dsolve({Eqns, ICs}, numeric); plots:-odeplot(soln, [r, fnl(r)], r0 .. S, color = col) end proc:

plotsol1s(-3*10^5, -1, 10^(-10), Rm, red)

Error, (in f) unable to store '[HFloat(0.0)]' when datatype=float[8]

 

``

``


Any help would be greatly appreciated.

 

Gambia Man

 

Hello

I need to solve the Bending Vibration of Euler-Bernouli Beam Problem and I keep getting stuck. I start with a fairly straight forward fourth order differential eq. Using the dsolve command gives me the general solution

Y(x)=A*sin(a*x)+B*cos(a*x)+C*sinh(a*x)+D*cosh(a*x)

Maple insist on using e^(x)+e^(-x) instead on sinh and cosh - but it's the same. So far so good.

My specific problem is a clamped-pinned beam of length l - so my boundary conditions are (correct me if I'm wrong here):

In the clamped end at x=0: Y(0)=0, Y'(0)=0

In the pinned end at x=l: Y(L)=0, Y''(0)=0

Using both the dsolve(ode,ics) and a dsolve(ode) and then solve(ics) both results in the trivial solution Y(x)=0 - which is wrong - I know there is a tan(a*l)-tanh(a*l) solution.

To get a easier and well documented example to solve by hand, I also tried with a simply supported beam. Boundary conditions are then:

Y(0)=Y(l)=0

Y''(0)=Y''(l)=0

Same result - only the trivial solution Y(x)=0 and If you solve it by hand you get a sin(a*l) solution.

 

What am I doing wrong? Is it syntax error on my part or what?

 

I have attached both my maple doc and a pdf with a walkthrough of the correct solution.

Beam_vibration.mw Transverse_vibration_of_beams.pdf

 

Any help would be appreciated

Kind regards

Jacob

 

and plot  function I? This I is the area which I wrote at the paper.

Could you give me the code which can be used to solve the ODE by numerical method and plot I with respect to t?

I think I have write down everything clearly but if you feel confused please ask me.

I am eager to know the code. Thanks very much!

Perhaps the question is trivial, but I could not find the solution.

I am solving numerically an ODE (e.g., a simple harmonic oscillator) with the righthand side that contains a random part. For example it is

eq:=diff(a(t),t,t) + a(t) = (1+0.01*R(t))*cos(5*t);

where R(t) is a random function of t. How can I make such a function?

The naive attempt: 

eq:=diff(a(t),t,t) + a(t) = 3.*(1+0.1*rand()/1000000000000.)*cos(5.*t);


gave me a fixed (while random) value, e.g.

...

But, I need this coefficeint to be random each time step.

Any suggestions are very welcome!

 

 

 

I want to plot E with respect to t.

I typed the sys of this equation and try to plot it with DEplot but failed.

This is the code

DEplot([sys], [x(t), y(t), E(t)], t = -25 .. 25, [[x(0) = -.3, y(0) = .2]], scene = [E(t), t], stepsize = 0.5e-1, linecolor = red, arrows = none);

The error is Error, (in DEtools/DEplot/CheckDE) derivatives must be given explicitly.

Eager to know the answer, please help me . Thanks very much.

 

Dear collagues

Hi,

I write a code to solve a system of ODE. It solve the ODES in a wide range of parameters but as I decrease NBT below 0.5, it doesnt converge. I do my best but I couldn't find the answer. Would you please help me? Thank you

Here is my code and it should be run for 0.1<NBT<10. the value of NBT is input directly in res1.

restart:
EPSILONE:=1000:
Digits:=15:

a[mu1]:=5.45:
b[mu1]:=108.2:
a[k1]:=1.292:
b[k1]:=-11.99:




rhop:=4175:
rhobf:=998.2:
mu1[bf]:=9.93/10000:
k1[bf]:=0.597:

rhost(eta):=1-phi(eta)+phi(eta)*rhop/rhobf;
k:=unapply(k1[bf]*(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2),eta);


eq1:=(diff(u(eta), eta))*a[mu1]*(diff(phi(eta), eta))+2*(diff(u(eta), eta))*b[mu1]*phi(eta)*(diff(phi(eta), eta))+((diff(u(eta), eta))+(diff(u(eta), eta))*a[mu1]*phi(eta)+(diff(u(eta), eta))*b[mu1]*phi(eta)^2)/(eta+EPSILONE)+diff(u(eta), eta, eta)+(diff(u(eta), eta, eta))*a[mu1]*phi(eta)+(diff(u(eta), eta, eta))*b[mu1]*phi(eta)^2+1-phi(eta)+phi(eta)*rhop/rhobf:
eq2:=(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2)*(diff(T(eta), eta))/(eta+EPSILONE) + (a[k1]*(diff(phi(eta), eta))+2*b[k1]*phi(eta)*(diff(phi(eta), eta)))*(diff(T(eta), eta))+(1+a[k1]*phi(eta)+b[k1]*phi(eta)^2)*(diff(T(eta), eta, eta));
eq3:=diff(phi(eta),eta)+phi(eta)/(N_bt*(1+gama1*T(eta))^2)*diff(T(eta),eta):

eq2:=subs(phi(0)=phi0,eq2):
eq3:=subs(phi(0)=phi0,eq3):


eval(dsolve({eq3,phi(1)=phiv},phi(eta)),(T)(1)=1):
Phi:=normal(combine(%)):
Teq:=isolate(eval(eq2,Phi),diff(T(eta),eta)):
ueq1:=eval(eq1,Phi)=0:
ueq2:=subs(Teq,ueq1):


lambda:=0;
Ha:=0;
N_bt:=cc*NBT+(1-cc)*0.8;
kratio:=k1[p]/k1[bf]:






GUESS:=[T(eta) =0.0001*eta, u(eta) =0.1*eta, phi(eta) = 0.3*(eta-1)^4];
res1 := dsolve(subs(NBT=0.48,gama1=0.2,phiv=0.06,{eq1,eq2,eq3,u(0)=lambda*D(u)(0),D(u)(1)=0,T(0)=0,phi(1)=phiv,T(1)=1}), numeric,method=bvp[midrich],maxmesh=4000,approxsoln=GUESS, output=listprocedure,continuation=cc):
G0,G1,G2:=op(subs(subs(res1),[phi(eta),u(eta),diff(T(eta),eta)])):

masst:=evalf(int((1-G0(eta)+G0(eta)*rhop/rhobf)*G1(eta), eta = 0..1));
heatt:=(1+a[k1]*G0(0)+b[k1]*G0(0)^2)*G2(0);

plots:-odeplot(res1,[[eta,T(eta)]],0..1,legend=[T],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);
plots:-odeplot(res1,[[eta,u(eta)]],0..1,legend=[u],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);
plots:-odeplot(res1,[[eta,phi(eta)]],0..1,legend=[phi],color=["Black"],linestyle=Solid,axes=boxed,thickness=3);

>
>

 

Thank you

 

Amir

Hello altogether,

I want to plot the numerical result of an ODE, which seems to be pretty simple at first sight, but the difficulty is that the boundaries are depending on the solution.

The following pseudo-code describes what I want to have, but it doesn't work. This code fills the RAM pretty fast and you will have to kill the process.

Free-boundary_prob.mw

Is it possible to calculate a solution to this problem numerically (or even analytically) and if yes, how?

Since I am new here, I am sorry for any bad-to-read maple code or any noob errors I have made. I would be very thankful, for any response and help.

Greetings

butterflyfart

EF.3.mwHi, I want to ask that how to find the exact solution of equation without applying any technique

Dear Collegues

I have a system of odes as follows

restart:
#gama1:=0.2:

#rhop:=5180:
#rhobf:=998.2:
#a[mu1]:=39.11:
#b[mu1]:=533.9:
#a[k1]:=7.47:
#b[k1]:=0:

Teq := N_bt*T(eta)^2*(exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta)))))^2*(diff(T(eta), eta, eta))*gama1^2*b[k1]+N_bt*T(eta)^2*exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta))))*(diff(T(eta), eta, eta))*gama1^2*a[k1]+2*N_bt*T(eta)*(exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta)))))^2*(diff(T(eta), eta, eta))*gama1*b[k1]+N_bt*T(eta)^2*(diff(T(eta), eta, eta))*gama1^2;


UEQ:=(a[mu1]*(-(diff(T(eta), eta))/(N_bt*(1+gama1)*(1+gama1*T(eta)))+(T(eta)-1)*gama1*(diff(T(eta), eta))/(N_bt*(1+gama1)*(1+gama1*T(eta))^2))*exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta))))+2*b[mu1]*(exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta)))))^2*(-(diff(T(eta), eta))/(N_bt*(1+gama1)*(1+gama1*T(eta)))+(T(eta)-1)*gama1*(diff(T(eta), eta))/(N_bt*(1+gama1)*(1+gama1*T(eta))^2)))*(diff(u(eta), eta))+(1+a[mu1]*exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta))))+b[mu1]*(exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta)))))^2)*(diff(u(eta), eta, eta))+1-exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta))))+exp(-(T(eta)-1)/(N_bt*(1+gama1)*(1+gama1*T(eta))))*rhop/rhobf;

I want to solve them with the following boundary conditions

T(0)=0, T(1)=1

u(0)=L*D(u)(0), D(u)(1)=0

However I tried, I cannot find the solution in a closed form. I want to know that is there a closed form solution for the above odes?

Thank you

Amir

When solving ODE with dsolve, I don't want to use "_C1" as the name of constant, I want to specify by myself, how can I do?

Hi,

In line with my previous questions

http://www.mapleprimes.com/questions/201944-Convergence-Problem-In-My-Algorithm-DSOLVE

I write a code with different variables and odes. However I tried to get the solution, I cannot get the results. I think it should be some problems with my Guess procedure. The code is attached. I would be most grateful if you help me on this problem.

code.mw

Thank you.

Amir

5 6 7 8 9 10 11 Last Page 7 of 30