List manipulation – Convert a backward / forward sweep code in Mathematica

I am trying to convert a simple backward / forward sweep code (with RK4 integrator) to Mathematica, but I get some error messages and cannot determine where I am coding incorrectly. The original code is written in MatLab and is in the reference "Maia Martcheva – An Introduction to Mathematical Epidemiology"::

function ocmodel1
% This function computes the optimal control
% and the corresponding solution using forward-backward ...
clear all;

test = -1;

 Δ = 0.001; %set tolerance
 N = 100; %number of subdivisions
 h = 1/N; %step
 t = 0:h:1; % t-variable mesh

 u = zeros(1,length(t)); %initialization
 x = zeros(1,length(t));
 lam = zeros(1,length(t));

 x(1) = 10; %initial value assigned to x(0)

 beta = 0.05; %parameters
 mu = 0.01;
 gamma = 0.5;
 P = 100;
 w1 = 1;

 while (test<0) % while the tolerance is reached, repeat
 oldu = u;
 oldx = x;
 oldlam = lam;

for i=1:N %loop that solve the forward ...
differential equation
k1 = beta*(P-x(i))*x(i) -(mu + gamma)*x(i) - ...
k2 = beta*(P-x(i)-0.5*k1*h)*(x(i)+0.5*k1*h) - ...
k3 = beta*(P-x(i)-0.5*k2*h)*(x(i)+0.5*k2*h) - ...
k4 = beta*(P-x(i)-k3*h)*(x(i)+k3*h) - ...

x(i+1) = x(i) + (h/6)*(k1+2*k2+2*k3+k4);


for i=1:N %loop that solves the backward ...
differential equation of the adjoint system
j = N + 2 -i;
k1 = ...
-w1-lam(j)*(beta*(P-x(j))-beta*x(j)-(mu+gamma) ...
- u(j));
k2 = ...
-w1-(lam(j)-0.5*k1*h)*(beta*(P-x(j)+0.5*k1*h) ...
-(mu+gamma) -0.5*(u(j)+u(j-1)));
k3 = ...
-w1-(lam(j)-0.5*k2*h)*(beta*(P-x(j)+0.5*k2*h) ...
-(mu+gamma) -0.5*(u(j)+u(j-1)));
k4 = -w1 -(lam(j)-k3*h)*(beta*(P-x(j)+k3*h) ...
-(mu+gamma) - u(j-1));

lam(j-1) = lam(j) - (h/6)*(k1+2*k2+2*k3+k4);


u1 = min(100,max(0,lam.*x/2));
u = 0.5*(u1 + oldu);

temp1 = Δ *sum(abs(u)) - sum(abs(oldu - u));
temp2 = Δ *sum(abs(x)) - sum(abs(oldx - x));
temp3 = Δ *sum(abs(lam)) - sum(abs(oldlam -lam));

test = min(temp1,min(temp2,temp3));


figure(1) %plotting



My attempt to write this code in Mathematica is as follows:

(* This function computes the optimal control
 % and the corresponding solution using forward-backward...
sweep *)
Clear( all)
test = -1;
Δ = 0.001;
n = 100;
h = 1/n;
t = Range(0, 1, h);
u = {};
x = {};
lam = {};
For(i = 1, i < n,
 AppendTo(u, 0);
 AppendTo(x, 0);
 AppendTo(lam, 0);
x = ReplacePart(x, 1 -> 10); (* initial value assigned to x(0) *)
beta = 0.05;(* parameters *)
mu = 0.01;
gamma = 0.5;
P = 100;
w1 = 1;
While(test < 0, (*  while the tolerance is reached,repeat*)
    oldu = u;
    oldx = x;
    oldlam = lam;
    For( i = 1, i < n, 
  k1 = beta*(P - x((i)))*x((i)) - (mu + gamma)*x((i)) - u((i))*x((i));
  k2 = (beta*(P - x((i))) - 0.5*k1*h)*(x((i)) + 0.5*k1*h) - (mu + 
       gamma)*(x((i)) + 0.5*k1*h) - 
    0.5*(u((i)) + u((i + 1)))*(x((i)) + 0.5*k1*h); 
  k3 = beta*(P - x((i)) - 0.5*k2*h)*(x((i)) + 0.5*k2*h) - (mu + 
       gamma)*(x((i)) + 0.5*k2*h) - 
    0.5*(u((i)) + u((i + 1)))*(x((i)) + 0.5*k2*h);
  k4 = beta*(P - x((i)) - k3*h)*(x((i)) + k3*h) - (mu + 
       gamma)*(x((i)) + k3*h) - u((i + 1))*(x((i)) + k3*h);
  ReplacePart(x, i + 1 -> x((i)) + (h/6)*(k1 + 2*k2 + 2*k3 + k4)); 
  i++ );
 For(i = 1, i < n, j = n + 2 - i; 
  k1 = -w1 - 
    lam((j))*(beta*(P - x((j))) - beta*x((j)) - (mu + gamma) - 
  k2 = -w1 - (lam((j)) - 
       0.5*k1*h)*(beta*(P - x((j)) + 0.5*k1*h) - (mu + gamma) - 
       0.5*(u((j)) + u((j - 1)))); 
  k3 = -w1 - (lam((j)) - 
       0.5*k2*h)*(beta*(P - x((j)) + 0.5*k2*h) - (mu + gamma) - 
       0.5*(u((j)) + u((j - 1)))); 
  k4 = -w1 - (lam((j)) - 
       k3*h)*(beta*(P - x((j)) + k3*h) - (mu + gamma) - u((j - 1)));
   j - 1 -> lam((j)) - (h/6)*(k1 + 2*k2 + 2*k3 + k4)); i++);
 u1 = Min(100, Max(0, lam.x/2));
 u = 0.5*(u1 + oldu);
 temp1 = Δ*Sum(Abs(u), {1, Length(u)}) - 
   Sum(Abs(oldu - u), {1, Length(u)});
 temp2 = Δ*Sum(Abs(x), {1, Length(x)}) - 
   Sum(Abs(oldx - x), {1, Length(x)});
 temp3 = Δ*Sum(Abs(lam), {1, Length(lam)}) - 
   Sum(Abs(oldlam - lam), {1, Length(lam)});

 test = Min(temp1, Min(temp2, temp3)); i++)
ListPlot(t, u)
ListPlot(t, x)

I get error messages about the indexes, but I can't see what I did wrong.

How do I run a Mathematica package (.m) from Python with the "Wolfram Client Library for Python"?

This question is about calling a Mathematica package (.m file) from Python using the "Wolfram Client Library for Python" (

For completeness, consider the following Mathematica package:

AddTwo::usage = "AddTwo(a, b) returns a+b";
AddTwo(a_, b_) := a + b;

target: Call Python's AddTwo from the Wolfram Client Library for Python. If this is not currently supported, add this feature to the library.

Documentation for the library can be found at:

A similar question asked before the Wolfram Client Library for Python was released is:
How do I run a Mathematica package (.m) from Python?

mathematica – Error installing VPN on Linux

I have an academic license to use mathematicaTherefore I have to be connected to my institute network via the VPN. I get the error during installation.
I tried to execute the following command: sudo apt-get install openvpn

Reading package lists... Done
Building dependency tree       
Reading state information... Done
openvpn is already the newest version (2.3.10-1ubuntu2.2).
The following packages were automatically installed and are no longer required:
  libdbusmenu-gtk4 linux-headers-4.4.0-128 linux-headers-4.4.0-128-generic
  linux-image-4.4.0-128-generic linux-image-extra-4.4.0-128-generic
Use 'sudo apt autoremove' to remove them.
0 upgraded, 0 newly installed, 0 to remove and 0 not upgraded.
1 not fully installed or removed.
After this operation, 0 B of additional disk space will be used.
Do you want to continue? (Y/n) y
Setting up matlab-support (0.0.21) ...
No default Matlab path found. Exiting.
dpkg: error processing package matlab-support (--configure):
 subprocess installed post-installation script returned error exit status 1
Errors were encountered while processing:
E: Sub-process /usr/bin/dpkg returned an error code (1)

A window always opens in which I have to define the path for matlab. Please help me with this

Fitting – NonlinearModelFit – Mathematica Stack Exchange

I tried to fit a model to the following experimental data:

data1 = {{0.00000000001, 100}, {1, 0.7}, {4, 0.03}, {24, 0.001}}

The first value should be 0.

NonlinearModelFit[data1, a x^-b, {a, b}, x]

gave the following equation: 0.390336/x^0.21896

The problem is that it doesn't match the data.
Using the same data set in Excel, the following equation was obtained: 0.62/x^-2.03
and this equation fits the data.
What is wrong with the Mathematica model fitting?

Calculus and Analysis – Why can't Mathematica return the answer for this integrand?

When I run Integrate Command for an expression does not return anything, as you can see in the screenshot below. What is the reason?
Enter the image description here

My input:

 1/(2 Sqrt(sp) (Beta)) E^(-re^2 (se - (Beta)^2/(4 sp)))
   Sqrt((Pi)) (-2 + 
    Erf((2 (re + rep) sp - re (Beta))/(2 Sqrt(sp))) + 
    Erf((-2 re sp + 2 rep sp + re (Beta))/(2 Sqrt(sp))) + 
    Erfc((2 (re + rep) sp + re (Beta))/(2 Sqrt(sp))) + 
    Erfc((2 rep sp - re (2 sp + (Beta)))/(2 Sqrt(sp)))), {re, 
  0, (Infinity)})

Fourth order tensor rotation – Mathematica Stack Exchange

Thank you for your reply to Mathematica Stack Exchange!

  • Please be sure answer the question. Provide details and share your research!

But avoid

  • Ask for help, clarify, or respond to other answers.
  • Make statements based on opinions; Support them with references or personal experiences.

Use MathJax to format equations. MathJax reference.

For more information, see our tips on writing great answers.

How does Mathematica calculate the gradient in the spherical coordinate system?

The case of a rectangular coordinate system

I already know that the differential equilibrium equations in different coordinate systems of elasticity can be calculated using a formula $ div ( sigma) + F = 0 $.

pdConv(f_) := 
  f /. Derivative(inds__)(g_)(vars__) :> 
    Apply(Defer(D(g(vars), ##)) &, 
     Transpose({{vars}, {inds}}) /. {{var_, 0} :> 
        Sequence(), {var_, 1} :> {var}}))
Thread(Div({{σ11(x, y, z), σ12(x, y, z), σ13(x, 
        y, z)}, {σ12(x, y, z), σ22(x, y, z), σ23(
        x, y, z)}, {σ13(x, y, z), σ23(x, y, 
        z), σ33(x, y, z)}}, {x, y, z}) + {F1(x, y, z), 
     F2(x, y, z), F3(x, y, z)} == 0) // pdConv

The physical equations in the rectangular coordinate system are as follows:

{{Subscript(ε, x), Subscript(ε, xy), 
    Subscript(ε, xz)}, {Subscript(ε, xy), 
    Subscript(ε, y), Subscript(ε, 
    yz)}, {Subscript(ε, xz), Subscript(ε, 
    yz), Subscript(ε, z)}} == 
   Grad({u(x, y, z), v(x, y, z), w(x, y, z)}, {x, y, z})) // Normal

The strain compatibility equation expressed by stress is as follows. Because of the symmetry, there are only six independent equations.

Curl(#, {x, y, z}) & /@ 
   Curl(#, {x, y, z}) & /@ {{e11(x, y, z), e12(x, y, z), e13(x, y, z)}, {e12(x, y, z), e22(x, y, z), e23(x, y, z)}, {e13(x, y, z), e23(x, y, z), e33(x, y, z)}}) // pdConv

The case of the spherical coordinate system

The physical equations in spherical coordinates are as follows:

{{εrr, εrθ, εrφ}, {εrθ, εθθ, εθφ}, {εrφ, εθφ, εφφ}} == 
    Grad({ur(r, θ, ϕ), uθ(r, θ, ϕ), 
      uφ(r, θ, ϕ)}, {r, θ, ϕ}, 
     "Spherical")) // Normal // pdConv

The differential equilibrium equations in spherical coordinates are as follows:

Thread(Div({{σrr(r, θ, φ), σrθ(r, θ, φ), σrφ(r, θ, φ)}, {σrθ(r, θ, φ), σθθ(r, θ, φ), σθφ(
        r, θ, φ)}, {σrφ( r, θ, φ), σθφ(r, θ, φ), σφφ(r, θ, φ)}}, {r, θ, φ}) + {Fr(r, θ, φ), Fθ(r, θ, φ), Fφ(r, θ, φ)} == 0) // pdConv

Problems to be solved

I want to know the specific algorithm of the physical equation under the specific spherical coordinate system (ie how to calculate Grad({ur(r, θ, ϕ), uθ(r, θ, ϕ), uφ(r, θ, ϕ)}, {r, θ, ϕ}, "Spherical")) specifically).

I also want to know how to solve the deformation compatibility equation in the spherical coordinate system with Mathematica.

Units – simplify unit rules? – Mathematica Stack Exchange

Code as follows:

In(1):= constants = Quantity /@ {"PlanckConstant", "BoltzmannConstant", "MolarGasConstant"};
  UnitConvert(Quantity("AvogadroNumber")*constants((2)))) // N);
  UnitSimplify /@ UnitConvert /@ constants // N


Out(3)= {Quantity(6.62607*10^-34, "Joules" "Seconds"), Quantity(1.38065*10^-23, ("Joules")/("Kelvins")), Quantity(8.31446, ("Kilograms" ("Meters")^2)/("Kelvins" "Moles" ("Seconds")^2))}

It seems that the unity of the molar gas constant cannot be simplified well enough, and very likely it was because of the "Moles" in unity when I calculate that molar gas constant by multiplying the Boltzmann constant With Avogadro's number Instead, the device can be simplified. There is no "Moles" in the unity of the latter.

In such cases, is there a way to improve automatic unit simplification? TKX.