restart;

with(plots): with(LinearAlgebra):

 

# TFSB Coefficients (symbolic in u)

beta0 := u -> (sin(u)*u^3 - 12*u^2 - 24*cos(u) + 24)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):

beta1 := u -> (5*sin(u)*u^3 + 12*cos(u)*u^2 + 24*cos(u) - 24)/(6*(sin(u)*u + 2*cos(u) - 2)*u^2):

beta2 := u -> beta0(u):

rho0 := u -> ((-u^2-12)*cos(u) - 5*u^2 + 12)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):

rho1 := u -> (-7*cos(u)*u^3 + 27*sin(u)*u^2 + 120*sin(u) - 120*u)/(60*u^2*(cos(u)*u + 2*u - 3*sin(u))):

rho2 := u -> -rho0(u):

 

# Secondary coefficients (simplified versions)

beta00 := u -> 13/42 - 9*u^2/7840:

beta10 := u -> 1/6 + u^2/720:

beta20 := u -> 1/42 - 17*u^2/70560:

beta01 := u -> 187/1680 + 611*u^2/705600:

beta11 := u -> 11/30 - 29*u^2/25200:

beta21 := u -> 37/1680 + 67*u^2/235200:

beta02 := u -> 11/70 + 491*u^2/352800:

beta12 := u -> 9/10 - 31*u^2/8400:

beta22 := u -> 31/70 + 811*u^2/352800:

 

rho01 := u -> 2/105 + 407*u^2/1058400:

rho11 := u -> -19/210 + 41*u^2/105840:

rho21 := u -> -1/168 - 101*u^2/529200:

rho02 := u -> 53/1680 + 1633*u^2/2116800:

rho12 := u -> 8/105 - 4*u^2/6615:

rho22 := u -> -101/1680 - 2273*u^2/2116800:

 

# Problem definition

omega := 1:

epsilon := 3*Pi/2:

phi := x -> 3*sin(x) - 5*cos(x):  # history function

 

f := (x, v, vp, vd) -> -v - vd + 3*cos(x) + 5*sin(x):

g := proc(x, v, vp, vd, vdp)

    local fx, fv, fvp, fvd;

    fx := -3*sin(x) + 5*cos(x);

    fv := -1;

    fvp := 0;

    fvd := -1;

    return fx + fv*vp + fvp*0 + fvd*vdp;

end proc:

 

# Initial conditions

a := 0: b := 10:

v0 := -5: vp0 := 3:

 

# Variable step-size parameters

tol := 1e-10:

h_min := 0.01:

h_max := 0.5:

h_init := Pi/8:

 

# Store results

X := [a]: V := [v0]: Vp := [vp0]:

h_curr := h_init:

x_curr := a:

v_curr := v0:

vp_curr := vp0:

 

# For history: need v at x-epsilon

get_v_delayed := proc(xx)

    if xx < a then return phi(xx);

    else

        # Interpolate from stored solution

        idx := 1;

        while idx < nops(X) and X[idx] < xx do idx := idx+1; end do;

        if idx = 1 then return phi(xx);

        elif X[idx] = xx then return V[idx];

        else

            # Linear interpolation

            return V[idx-1] + (V[idx]-V[idx-1])*(xx-X[idx-1])/(X[idx]-X[idx-1]);

        end if;

    end if;

end proc:

 

# Newton solver for block

solve_block := proc(x0, v0, vp0, h, omega)

    local u, bet0, bet1, bet2, rho0, rho1, rho2, bet00, bet10, bet20, bet01, bet11, bet21, bet02, bet12, bet22,

          rho01, rho11, rho21, rho02, rho12, rho22, F, J, V0, V1, V2, Vp0, Vp1, Vp2, tolN, iter, dv, dV;

   

    u := omega*h;

    bet0 := beta0(u); bet1 := beta1(u); bet2 := beta2(u);

    rho0 := rho0(u); rho1 := rho1(u); rho2 := rho2(u);

    bet00 := beta00(u); bet10 := beta10(u); bet20 := beta20(u);

    bet01 := beta01(u); bet11 := beta11(u); bet21 := beta21(u);

    bet02 := beta02(u); bet12 := beta12(u); bet22 := beta22(u);

    rho01 := rho01(u); rho11 := rho11(u); rho21 := rho21(u);

    rho02 := rho02(u); rho12 := rho12(u); rho22 := rho22(u);

   

    # Initial guesses

    V1 := v0 + h*vp0;

    V2 := v0 + 2*h*vp0;

    Vp1 := vp0;

    Vp2 := vp0;

   

    tolN := 1e-12;

    for iter from 1 to 10 do

        # Compute delayed values

        vd0 := get_v_delayed(x0 - epsilon);

        vd1 := get_v_delayed(x0 + h - epsilon);

        vd2 := get_v_delayed(x0 + 2*h - epsilon);

        vdp0 := (get_v_delayed(x0 - epsilon + 1e-8) - vd0)/1e-8;

        vdp1 := (get_v_delayed(x0 + h - epsilon + 1e-8) - vd1)/1e-8;

        vdp2 := (get_v_delayed(x0 + 2*h - epsilon + 1e-8) - vd2)/1e-8;

       

        # Compute gamma and g

        gam0 := f(x0, v0, vp0, vd0);

        gam1 := f(x0+h, V1, Vp1, vd1);

        gam2 := f(x0+2*h, V2, Vp2, vd2);

        g0 := g(x0, v0, vp0, vd0, vdp0);

        g1 := g(x0+h, V1, Vp1, vd1, vdp1);

        g2 := g(x0+2*h, V2, Vp2, vd2, vdp2);

       

        # Residuals

        F1 := h*vp0 - (V1 - v0 + h^2*(bet00*gam0 + bet10*gam1 + bet20*gam2)

              + h^3*(rho01*g0 + rho11*g1 + rho21*g2));

        F2 := h*Vp1 - (V1 - v0 + h^2*(bet01*gam0 + bet11*gam1 + bet21*gam2)

              + h^3*(rho01*g0 + rho11*g1 + rho21*g2));

        F3 := h*Vp2 - (V1 - v0 + h^2*(bet02*gam0 + bet12*gam1 + bet22*gam2)

              + h^3*(rho02*g0 + rho12*g1 + rho22*g2));

        F4 := V2 - (2*V1 - v0 + h^2*(bet0*gam0 + bet1*gam1 + bet2*gam2)

              + h^3*(rho0*g0 + rho1*g1 + rho2*g2));

       

        F := Vector([F1, F2, F3, F4]);

        if LinearAlgebra:-Norm(F) < tolN then break; end if;

       

        # Approximate Jacobian (finite differences)

        J := Matrix(4,4);

        delta := 1e-6;

        for j from 1 to 4 do

            V_pert := Vector([V1, V2, Vp1, Vp2]);

            V_pert[j] := V_pert[j] + delta;

            V1p := V_pert[1]; V2p := V_pert[2]; Vp1p := V_pert[3]; Vp2p := V_pert[4];

            gam1p := f(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon));

            gam2p := f(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon));

            g1p := g(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon),

                     (get_v_delayed(x0+h-epsilon+1e-8)-get_v_delayed(x0+h-epsilon))/1e-8);

            g2p := g(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon),

                     (get_v_delayed(x0+2*h-epsilon+1e-8)-get_v_delayed(x0+2*h-epsilon))/1e-8);

           

            F1p := h*vp0 - (V1p - v0 + h^2*(bet00*gam0 + bet10*gam1p + bet20*gam2p)

                   + h^3*(rho01*g0 + rho11*g1p + rho21*g2p));

            F2p := h*Vp1p - (V1p - v0 + h^2*(bet01*gam0 + bet11*gam1p + bet21*gam2p)

                   + h^3*(rho01*g0 + rho11*g1p + rho21*g2p));

            F3p := h*Vp2p - (V1p - v0 + h^2*(bet02*gam0 + bet12*gam1p + bet22*gam2p)

                   + h^3*(rho02*g0 + rho12*g1p + rho22*g2p));

            F4p := V2p - (2*V1p - v0 + h^2*(bet0*gam0 + bet1*gam1p + bet2*gam2p)

                   + h^3*(rho0*g0 + rho1*g1p + rho2*g2p));

           

            Fp := Vector([F1p, F2p, F3p, F4p]);

            J[1..4, j] := (Fp - F)/delta;

        end do;

       

        dV := LinearAlgebra:-LinearSolve(J, -F);

        V1 := V1 + dV[1]; V2 := V2 + dV[2]; Vp1 := Vp1 + dV[3]; Vp2 := Vp2 + dV[4];

    end do;

   

    return [V1, V2, Vp1, Vp2];

end proc:

 

# Main variable step-size loop

printf("Variable step-size integration for Example 1\n");

printf("tol = %e, h_init = %f\n", tol, h_init);

 

while x_curr < b - 1e-12 do

    # Try current step

    sol := solve_block(x_curr, v_curr, vp_curr, h_curr, omega);

    V1 := sol[1]; V2 := sol[2]; Vp1 := sol[3]; Vp2 := sol[4];

   

    # Compute with two half-steps

    sol_half1 := solve_block(x_curr, v_curr, vp_curr, h_curr/2, omega);

    V_mid := sol_half1[2]; Vp_mid := sol_half1[4];

    sol_half2 := solve_block(x_curr + h_curr/2, V_mid, Vp_mid, h_curr/2, omega);

    V2_half := sol_half2[2];

   

    # Error estimate

    err := abs(V2 - V2_half) / (2^6 - 1);

   

    if err < tol then

        # Accept step

        x_next := x_curr + 2*h_curr;

        X := [op(X), x_curr + h_curr, x_next];

        V := [op(V), V1, V2];

        Vp := [op(Vp), Vp1, Vp2];

        x_curr := x_next;

        v_curr := V2;

        vp_curr := Vp2;

        printf("x = %7.4f, h = %8.5f, err = %12.5e\n", x_curr, h_curr, err);

       

        # Adjust step size

        if err < tol/2 then

            h_curr := min(2*h_curr, h_max);

        end if;

    else

        # Reject step, reduce h

        h_curr := max(h_curr/2, h_min);

        printf("  Rejecting, new h = %8.5f\n", h_curr);

    end if;

end do:

 

# Exact solution for comparison

exact := x -> 3*sin(x) - 5*cos(x);

errors := [seq(abs(V[i] - exact(X[i])), i=1..nops(X))];

 

# Visualization

p1 := pointplot([seq([X[i], errors[i]], i=1..nops(X))], color=red, symbol=circle,

                title="Example 1: Variable Step-Size TFSB - Absolute Errors",

                labels=["x", "Error"], labeldirections=[horizontal,vertical]);

p2 := plot([[x_curr, h_curr]], x=a..b, style=point, color=blue,

            title="Step-size evolution", labels=["x", "h"]);

display(p1);

display(p2);

 

printf("\nFinal results for Example 1:\n");

printf("Number of steps: %d\n", nops(X)-1);

printf("Maximum error: %e\n", max(errors));

printf("Final step-size: %f\n", h_curr);

I have Maple 2016 and Matlab 2016b installed on my Windows destop.

When I run

with(Matlab);
openlink()

in a Maple worksheet the following error window appears.

After closing the window I also see the following error message,

Error, (in Matlab:-openlink) there was a problem finding or loading matlink.dll. Refer to ?Matlab,setup for help configuring your system to work with the Matlab-link.

How to fix the connection from Maple to Matlab?

Good day, everyone

I am trying to code HPM but it's giving me the error code "Error, invalid subscript selector
" once I increase the number of iterations above 2. 

Attached below is the code. 

HPM.mw

Thanks 

Good day, all.

Please, I am working on the following code but found out that the command solve is not displaying any result. Your assistance and suggestions would be appreciated. Thank you, and best regards.

 

restart;
NULL;
t := sum(a[j]*q^j, j = 0 .. 9);
H := diff(t, q);
F := diff(t, q $ 2);
p1 := simplify(eval(t, q = x)) = y[n];
p2 := simplify(eval(F, q = x)) = f[n];
p3 := simplify(eval(F, q = x + h/4)) = f[n + 1/4];
p4 := simplify(eval(F, q = x + h/2)) = f[n + 1/2];
p5 := simplify(eval(F, q = x + (3*h)/4)) = f[n + 3/4];
p6 := simplify(eval(F, q = x + h)) = f[n + 1];
p7 := simplify(eval(F, q = x + (5*h)/4)) = f[n + 5/4];
p8 := simplify(eval(F, q = x + (3*h)/2)) = f[n + 3/2];
p9 := simplify(eval(F, q = x + (7*h)/4)) = f[n + 7/4];
p10 := simplify(eval(F, q = x + 2*h)) = f[n + 2];
r := seq(a[i], i = 0 .. 9);
s := p || (1 .. 10);

solve({s}, {r});

i used: 
Y := ssystem("dir C:"); print(Y)
result is

[0,"\" El volumen de la unidad C es OS\n El n£mero de serie del volumen es: 54A9-09DA\n\n Directorio de C:\\Program Files\\Maple 2016"

Windows operating system English version but Maple shows Spanish result
I want result is English. Please help me

(1)

Download language_maple.mw

How can I display the symbol (blue solid circle) and the line (blue line) together in the legend box (( from   to ))?

>

Download Plot1.mw

I have a file TEST.m. How can I make it so that every time I start Maple, all the subprograms in the TEST.m file will run first? Then I just need to type the function with(TEST): sumpro(2,3,4) to get the result 9. I copied the TEST.m file into Maple's lib directory, but it doesn't run after starting Maple.

I just need to type sumvip(2, 3, 4) to get the result, but Maple doesn't understand it.

Please help.

(1) sum of 2 , 3 and 4 is 9

Download TEST.mw

Assuming I have a sumpro function written in Maple 2016. How can I implement it in C# and what is the process?

Please help me.

sumpro := proc (i) local a, b;

        a := (rand(1 .. 10))(); b := (rand(1 .. 10))();

        print("Sum of ", a, " and ", b, " is ", a+b)

end proc;

save sumpro, "D://Sumpro.m"

Download mapleprime_sumpro_to_c.mw

Hi

How we can reduce run time for Cases III and IV ? (Maple 2016)

Thanks alot

Formula_II.mw

Good day, all.

Please I want to solve the following delay differential equation:

ODE := diff(y(t), t$2) = (2*(1-y(t-1)^2))*(diff(y(t), t))-y(t)

ics := y(0) = 1, (D(y))(0) = 0

using the following codes but there is an error. Please kindly help to modify the codes.

restart;
Digits:=30:

f:=proc(n)
	2*(1-(y[n-1])^2)*delta[n]+y[n]:
end proc:

g:=proc(n)
	-4*y[n-1]*delta[n-1]+2*(1-(y[n-1])^2)*f(n)-delta[n]:
end proc:


e1:=y[n+2] = -y[n]+2*y[n+1]+(1/120)*h^2*(-3*h*g(n+2)+3*g(n)*h+16*f(n+2)+16*f(n)+88*f(n+1)):
e2:=h*delta[n] = -y[n]+y[n+1]-(1/1680)*h^2*(-128*h*g(n+1)-11*h*g(n+2)+59*g(n)*h+40*f(n+2)+520*f(n)+280*f(n+1)):
e3:=h*delta[n+1] = -y[n]+y[n+1]+(1/1680)*h^2*(-152*h*g(n+1)-10*h*g(n+2)+32*g(n)*h+37*f(n+2)+187*f(n)+616*f(n+1)):
e4:=h*delta[n+2] = -y[n]+y[n+1]+(1/1680)*h^2*(128*h*g(n+1)-101*h*g(n+2)+53*g(n)*h+744*f(n+2)+264*f(n)+1512*f(n+1)):

inx:=0:
ind:=0:
iny:=1:
h:=1/2:
n:=1:
omega:=10:
u:=omega*h:
N:=solve(h*p = 10, p):

err := Vector(round(N)):
exy_lst := Vector(round(N)):
numerical_y1:=Vector(round(N)):

c:=1:
for j from 0 to 2 do
	t[j]:=inx+j*h:
end do:

vars:=y[n+1],y[n+2],delta[n+1],delta[n+2]:

step := [seq](eval(x, x=c*h), c=1..N):
printf("%6s%45s%45s\n", 
	"h","Num.y","Num.z");
#eval(<vars>, solve({e||(1..4)},{vars}));


st := time():
for k from 1 to N/2 do

	par1:=x[0]=t[0],x[1]=t[1],x[2]=t[2]:
	par2:=y[n]=iny,delta[n]=ind:
	res:=eval(<vars>, fsolve(eval({e||(1..4)},[par1,par2]), {vars}));

	for i from 1 to 2 do
		printf("%6.5f%45.30f%45.30f\n", 
		h*c,res[i],res[i+2]):
		
		numerical_y1[c] := res[i]:
		
		c:=c+1:
	end do:
	iny:=res[2]:
	ind:=res[4]:
	inx:=t[2]:
	for j from 0 to 2 do
		t[j]:=inx + j*h:
	end do:
end do:
v:=time() - st;
v/4;
printf("Maximum error is %.13g\n", max(err));
NFE=evalf((N/4*3)+1);
#get array of numerical and exact solutions for y1
numerical_array_y1 := [seq(numerical_y1[i], i = 1 .. N)]:
#exact_array_y1 := [seq(exy[i], i = 1 .. N)]:

#get array of time steps
time_t := [seq](step[i], i = 1 .. N):

#display graphs for y1
with(plots):
numerical_plot_y1 := plot(time_t, numerical_array_y1, style = [point], symbol = [asterisk],
				color = [blue,blue],symbolsize = 20, legend = ["TFIBF"]);

 Thank you, and best regards.

How can numbers be displayed inside the contour plot?

 restart;
with(plots);
contourplot(x*exp(-x^2 - y^2), x = -2 .. 2, y = -2 .. 2, axes = boxed);
like this

I need admin's help
I use evalf(3*21/100,3)=0.630
and evalf(3*89/100,3)=2.67
Is there a way for me to get 2 decimal places
so evalf(3*21/100,3)=0.63?

Good day everyone,

I am writing a numerical code using dsolve which works fine but I have a challenge in inputting the previous answers in the subsequent ones. For example, how can I substitute the solutions in S1 into equ11, equ22, equ33, and equ44 in the link below? 

Thank you very much as I will be expecting responses from you soon.

New.mw

I want to express my two variable function f using Taylor expansion. But no success yet.

Why Taylor series can not estimate my function in desired interval [-1<x,y<1]?

>  >  with(Student[MultivariateCalculus]):   >  >  >  >    >  >  plot3d(h, xx = -1 .. 1, yy = -1 .. 1, color = red, style = surface) >

Download taylorProblem.mw

How will I use maple 2016 to solve ODEs and showing the steps involved because this will increase my understanding in it.