Hi Users!

I hope everyone is fine here. I have plotted the density plot below:

restart; Digits := 20; with(LinearAlgebra); with(plots); N := 20; Mx := 20; L := 1; `&Delta;x` := L/Mx; T := 2; `&Delta;t` := T/N; for i from 0 while i <= Mx do u[i, 0] := 0 end do; for n from 0 while n <= N do u[0, n] := 0; u[Mx, n] := 0 end do; for n from 0 while n <= N-1 do for i while i <= Mx-1 do Ru[i, n] := eval((u[i, n+1]-u[i, n])/`&Delta;t` = (u[i+1, n+1]-2*u[i, n+1]+u[i-1, n+1])/`&Delta;x`^2-u[i, n+1]+.5) end do; Sol[n] := fsolve({seq(Ru[i, n], i = 1 .. Mx-1)}); assign(op(Sol[n])) end do;

Digits := 10; NP := 100; XX := [seq(seq(i1/Mx, i1 = 0 .. Mx), i2 = 0 .. N)]; TT := [seq(seq(i2, i1 = 0 .. Mx), i2 = 0 .. N)]; ZZ := [seq(seq(u[i1, i2], i2 = 0 .. N), i1 = 0 .. Mx)]; interfunc := subs(__M = Matrix(Matrix(Mx+1, N+1, ZZ), datatype = float[8]), proc (x, y) options operator, arrow; CurveFitting:-ArrayInterpolation([[`$`(i1, i1 = 0 .. Mx)], [`$`(i2, i2 = 0 .. N)]], __M, [[x], [y]], method = cubic)[1, 1] end proc); newz := CurveFitting:-ArrayInterpolation([[`$`(i1, i1 = 0 .. Mx)], [`$`(i2, i2 = 0 .. N)]], Matrix(Mx+1, N+1, ZZ, datatype = float[8]), [[seq(Mx*(i3-1)/(NP-1), i3 = 1 .. NP)], [seq(N*(i3-1)/(NP-1), i3 = 1 .. NP)]], method = cubic); nminz, nmaxz := (min, max)(newz); C := .666*(1-ImageTools:-FitIntensity(newz)); PC := PLOT(GRID(0 .. Mx, 0 .. N, newz, COLOR(HUE, C)), STYLE(PATCHNOGRID)); numcontours := 15; PP := (proc (P) options operator, arrow; (op(0, P))(op(P), ROOT(BOUNDS_X(0), BOUNDS_Y(0), BOUNDS_WIDTH(600), BOUNDS_HEIGHT(500))) end proc)(plots:-display(PC, plots:-contourplot(interfunc, 0 .. Mx, 0 .. N, thickness = 0, contours = [seq(nminz+(nmaxz-nminz)*(i3-1)/(numcontours+2-1), i3 = 1 .. numcontours+2)]), seq(plot(ZZ[1], nminz .. nminz, thickness = 15, color = COLOR(HUE, .666*(1-i3/(numcontours+1))), legend = sprintf(" %.3f", nminz+i3*(nmaxz-nminz)/(numcontours+1))), i3 = numcontours+1 .. 0, -1), legendstyle = [location = right, font = [Helvetica, 14]], font = [Helvetica, 16], labelfont = [Helvetica, bold, 16], labels = [x, t], labeldirections = [horizontal, vertical]));
plots[display](PP, size = [500, 400]);

Here the t-axis is from 0 to 20 but its actual value is from 0 to 2 and the x-axis is from 0 to 20 but its actual value is from 0 to 2. How can I change the axis? Moreover, I used the following way to extract the data in a dat file to plot the function (say f) in some professional software.

with(plots); f := plot(sin(x), x = -Pi .. Pi); dat1 := `~`[plottools:-getdata]([f]); for i while i <= 1 do A[i] := dat1[i, 3]; Y[i] := A[i][1 .. -1, 2] end do; X := A[1][1 .. -1, 1]; MM1 := `<|>`(X, `$`(Y[j2], j2 = 1 .. 1)); ExportMatrix("C:/Users/Usman/Desktop//Graph f.dat", MM1);

How I can, I extract data in the form of a dat file to plot in some professional software? 

(I would prefer a solution for Maple 2015, but answers relative to newer versions are welcome)

Is there a simple way to force the result -y(1) + y(2) without using one of these two tricks?

# how can I get the expression of
int(diff(y(x), x), x=1..2);
                      / d                  \
                   int|--- y(x), x = 1 .. 2|
                      \ dx                 /

# Trick 1
int(diff(y(x), x), x);
eval(%, x=2)-eval(%, x=1)
                              y(x)
                          -y(1) + y(2)

# Trick 2
J := Int(diff(y(x), x), x = 1..2): 
value(IntegrationTools:-Parts(J, 1));
                          -y(1) + y(2)

TIA

Dear Users!

I hope everyone is fine. I want to plot the following sequence in 3d for t=0..1 and x=-pi..pi;

[0., 0.4995839572e-1*sin(x), 0.9966865249e-1*sin(x), .1488899476*sin(x), .1973955598*sin(x), .2449786631*sin(x), .2914567945*sin(x), .3366748194*sin(x), .3805063771*sin(x), .4228539261*sin(x), .4636476090*sin(x), .5028432109*sin(x), .5404195003*sin(x), .5763752206*sin(x), .6107259644*sin(x), .6435011088*sin(x), .6747409422*sin(x), .7044940642*sin(x), .7328151018*sin(x), .7597627549*sin(x), .7853981634*sin(x)];

In the sequence first entry (0) for t=0, second (0.4995839572e-1*sin(x)) for t=0.05, third (0.9966865249e-1*sin(x)) for t= 0.1 and so on the last entry (.7853981634*sin(x)) for t=1. In addition, how do I plot if the number of points exceeds in the sequence for example 100 or 1000 points, but the difference between two consecutive values for t is the same here the difference is Delta*t=0.05.

Note: I would -prefer an answer for Maple 2015, but I can accommodate an answer for a more recent version.

I have a function Gpdf from IR² to IR⁺ of class C¹ (this comes from the way this function is built).
Although its level curves are continuous, their display show discontinuities for some level values. 

The reason is that  Gpdf contains a term whose denominator vanishes and so, even if the left and right limits of Gpdf are the same at the vanishing point, the resulting plot is dicontinuous.

More details are given in the attached file Discontinuous_contours.mw.

I have tried to adjust the plotting grid, or even to superimpose contours drawn in domains containing no singularities, but I wasn't capable to get continuous drawings (see the attached file).

Do you have any idea to achieve this?

TIA

Hi!

I have M number of linear differential equations. I have solved this system using the 1,2,3,4 stag RK method in the attached file but did not find a significant difference in the accuracy. Kindly see what's wrong there. 

Thank you!

s-stage.mw

Dear Users!

I hope you are doing well. I have the following discretized form

for n>=1 and j=0..M. We obtained the following matrix equation for any "n" and j=0..M as:

I want matrix proc of any useful way to define A^n, u^n, and b^n. I am waiting for your positive response. Thanks in advancs

Dear Users!

I hope everyone is fine here. I have three vectors say V[1], V[2] and V[3] as:

restart; with(LinearAlgebra); with(linalg);
V[1] := Vector([3, 2, -4]);
V[2] := Vector([1, 2/3, 7]);
V[3] := Vector([-9, 0, 1/2]);

I can define a Vector V have V[1], V[2] and V[3] using blockmatrix command as:

C := blockmatrix(3, 1, [V[1], V[2], V[3]]);

My question is that how I can extract vectors V[1], V[2] and V[3] from C? Thanks in advance

Dear Users!

I hope you are doing well. In the attached file I want to convert the system of ODEs (attained the system against the value of M) into matrix form and need the matrices A, B and vector b. Remember the order of A, B and b vary as M vary. I am waiting for your kind response. Please take care and thanks

System_of_ODEs.mw

Dear Users!

I hope everyone is fine here. I want to solve the following system of PDEs associated with Robin-type boundary conditions. But got the error. Kindly help me to fix this issue. Thanks

restart; TT := 0.1e-2; l := 1/5; b[1] := .18; b[2] := 2*10^(-9); k[1] := 1.3*10^(-7); k[-1] := 24; k[2] := 7.2; p := .9997; d[1] := 0.412e-1; f := .2988*10^8; g := 2.02*10^7; s := 1.36*10^4; E[0] := 3.3*10^5; T1[0] := .5*10^9; C1[0] := 3.3*10^5; alpha[0] := 10^(-10); D1 := 10^(-6); D2 := 10^(-2); D3 := 10^(-6); d[4] := 1.155*10^(-2); t[0] := 1/D1; kappa := 10^4; k[3] := 300*(24*60); chi := 0; sigma := d[1]*t[0]; rho := f*t[0]*C1[0]/(E[0]*T1[0]); mu := k[1]*t[0]*T1[0]; eta := g/T1[0]; epsilon := t[0]*C1[0]*(p*k[2]+k[-1])/E[0]; omega := D3/D1; beta1 := b[1]*t[0]; beta2 := b[2]*T1[0]; phi := k[1]*t[0]*E[0]; lambda := t[0]*C1[0]*(k[-1]+k[2]*(1-p))/T1[0]; psi := t[0]*(k[-1]+k[2]); gamma1 := chi*alpha[0]/D1; delta := D2/D1; kappa := k[3]*t[0]*C1[0]/alpha[0]; xi := d[4]*t[0]; PDE1 := diff(u(y, t), t) = diff(u(y, t), y, y)-gamma1*(u(y, t)*(diff(theta(y, t), y, y))+(diff(u(y, t), y))*(diff(theta(y, t), y)))+sigma*piecewise(y <= l, 0, 1)+rho*C(y, t)/(eta+T(y, t))-sigma*u(y, t)-mu*u(y, t)*T(y, t)+epsilon*C(y, t); PDE2 := diff(theta(y, t), t) = delta*(diff(theta(y, t), y, y))+kappa*C(y, t)-xi*theta(y, t); PDE3 := diff(T(y, t), t) = omega*(diff(T(y, t), y, y))+beta1*(1-beta2*T(y, t))*T(y, t)-phi*u(y, t)*T(y, t)+lambda*C(y, t); PDE4 := diff(C(y, t), t) = mu*u(y, t)*T(y, t)-psi*C(y, t); ICs := u(y, 0) = piecewise(0 <= y and y <= l, 0, 1-exp(-1000*(x-l)^2)), T(y, 0) = piecewise(0 <= y and y <= l, 1-exp(-1000*(x-l)^2), 0), C(y, 0) = piecewise(l-epsilon <= y and y <= l+epsilon, exp(-1000*(x-l)^2), 1-exp(-1000*(x-l)^2)), theta(y, 0) = 0; BCs := {(D[1](C))(0, t) = 0, (D[1](C))(1, t) = 0, (D[1](T))(0, t) = 0, (D[1](T))(1, t) = 0, (D[1](theta))(0, t) = 0, (D[1](theta))(1, t) = 0, (D[1](u))(0, t) = 0, (D[1](u))(1, t) = 0};

PDE:= {PDE1, PDE2, PDE3, PDE4}; pds := pdsolve(PDE, {ICs}, BCs, numeric, spacestep = 1/100, timestep = 1/100);

Error, (in pdsolve/numeric/process_PDEs) specified dependent variable(s) {(D[1](C))(0, t) = 0, (D[1](C))(1, t) = 0, (D[1](T))(0, t) = 0, (D[1](T))(1, t) = 0, (D[1](theta))(0, t) = 0, (D[1](theta))(1, t) = 0, (D[1](u))(0, t) = 0, (D[1](u))(1, t) = 0} not present in input PDE
 

Hi!

I hope everyone is fine. I have a square matrix like the following form
A := Matrix([[10, -1, 2, 0], [-1, 11, -1, 3], [2, -1, 10, -1], [0, 3, -1, 8]]);
How to split A into three matrices D, L and U as:

D:= Matrix([[10, 0, 0, 0], [0, 11, 0, 0], [0, 0, 10, 0], [0, 0, 0, 8]]);
L := Matrix([[0, 0, 0, 0], [-1, 0, 0, 0], [2, -1, 0, 0], [0, 3, -1, 0]]);
U := Matrix([[0, -1, 2, 0], [0, 0, -1, 3], [0, 0, 0, -1], [0, 0, 0, 0]]);

I am waiting for your positive response. Please take care

Hi dear Users!

I hope everyone here is fine. I have a function like

f := exp(-t)*(x^2-5*x^3+10*x^5+x+3+.5*x^4)+(1/2)*x^2*(x-1)+2*sin(x);

I have to find the value of t at which the behavior of this function is constant for 6 decimal places against x from 0..1. This is my effort

restart;
N := 20; TOL := 10^(-6); Points := 100000;
f := exp(-t)*(x^2-5*x^3+10*x^5+x+3+.5*x^4)+(1/2)*x^2*(x-1)+2*sin(x);
for j from 0 while j <= 10 do print("\nWhen x = ", j/(10.));
for i from 0 while i <= Points do
g[i, j] := evalf(eval(f, [x = (1/10)*j, t = N*i/Points]));
if `and`(i >= 1, abs(g[i, j]-g[i-1, j]) < TOL) then print("Value of t = ", evalf(N*i/Points)); print("Value of f = ", g[i, j]); break else  
end if end do end do;

The same value is then verified by making graphs

plot([eval(f, t = 1), eval(f, t = 2), eval(f, t = 3), eval(f, t = 4), eval(f, t = 5), eval(f, t = 6), eval(f, t = 7)], x = 0 .. 1, color = [red, green, blue, cyan, yellow, black, purple]);
plot(eval(f, x = .8), t = 0.1e-1 .. N);

Here I want to know if is there any more effective maple command to find the value of t rather than using procedures (highlighted by red) or a graphical way.

Hi Users!

I hope everyone is fine. I want to plot any function say
f := exp(cos(x)+sin(x)) for x=a..b for any n say 12 so that h := (b-a)/n. 

For a=0, b=3 and n=12 I got h=1/4 and plot of f is:

But I want the plotting as given bellow where the value of f(x) is mentioned and girds line.

I am waiting for your answer. Thanks in advance.

Dear Users!

I hope everyone is fine here. I wrote the following statements with the print command:

restart; NN := [4, 6, 8]; a := 0; b := 2; n := 4;
h := evalf((b-a)/n); print("The domain of intergation is [a,b] = ", [a, b]);
f := exp(x); print("The given function is ", f);
Exact := evalf(int(f, x = a .. b)); print("The exact integration in [a,b] is ", Exact);
print("The value of h to divide the domain [a,b] into n subintervals is ", h);
print("Numerical integration in [a,b] is going to perform when h via RECTANGULAR METHOD for n = ", n);

The output is:

                          0.5000000000
        "The domain of integration is [a,b] = ", [0, 2]
                "The given function is ", exp(x)
       "The exact integration in [a,b] is ", 6.389056099
"The value of h to divide the domain [a,b] into n subintervals is ", 0.5000000000
"Numerical integration in [a,b] is going to perform when h via RECTANGULAR METHOD for n = ", 4

I want the actual values of a,b, n and h highlighted in above as:

                            0.5000000000
        "The domain of integration is [a,b] = ", [0, 2]
                "The given function is ", exp(x)
       "The exact integration in [0,2] is ", 6.389056099
"The value of h to divide the domain [0,2] into 4 subintervals is ", 0.5000000000
"Numerical integration in [0,2] is going to perform when 0.5 via RECTANGULAR METHOD for n = ", 4

Dear users! 

I hope everyone is fine here. I have the following expression:

r*y[0, 1]+y[0, 0]+(1/6)*r*(r-1)*(1+r)*y[-1, 3]+(1/2)*r*(r-1)*y[-1, 2]+(1/120)*r*(r-1)*(r-2)*(1+r)*(2+r)*y[-2, 5]+(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+(1/720)*r*(r-1)*(r-2)*(r-3)*(r-4)*(1+r)*y[-3, 8]+(1/360)*r*(r-1)^2*(r-2)*(r-3)*(1+r)*y[-3, 7]+(1/720)*r*(r-1)*(r-2)*(r-3)*(1+r)*(2+r)*y[-3, 6]:

What is the procedure to select some terms in the above expression for example for N=2 I just want the following terms:

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2];

and for N=3 I just want the following terms:

y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3];

and for N=4 I want:

(1/24)*r*(r-1)*(r-2)*(1+r)*y[-2, 4]+y[0, 0]+r*y[0, 1]+(1/2)*r*(r-1)*y[-1, 2]+(1/6)*r*(r-1)*(1+r)*y[-1, 3];

and so on,

in the descending order in first suffices and ascending order in second suffices (like term having y[0,0],  y[0,1], y[-1,2], y[-1,3], y[-2,4]). I am waiting for your response. Thanks.

How to force Maple to prove equality (2) under conditions cond.

Download Stable.mw

Thanks for your help.