0
$\begingroup$

I have a $N \times 1$ vector $F(t)$ for $0\leq t \leq 1$.

I want to find a $N\times N$ Matrix $C$ (a real-valued matrix) satisfying the following equation: $$\int_0^t F(t_0)dt_0= C.F(t)$$ If it is not possible,

how to find $C$ satisfying $\int_0^t F(t_0)dt_0 \simeq C.F(t)$ with minimum error?

For Example; suppose that we have a vector $F(t)$ as follows:

F := Function[{t}, {Piecewise[{{Sqrt[2], 0 <= t <= 1/2}}, 0], 
    Piecewise[{{2^(Sqrt[3]*(-(1/2) + 2*t)), 0 <= t <= 1/2}}, 0], 
       Piecewise[{{2^(3*Sqrt[5]*(1/6 - 2*t + 4*t^2)), 
       0 <= t <= 1/2}}, 0], Piecewise[{{Sqrt[2], 1/2 <= t <= 1}}, 0], 
       Piecewise[{{2^(Sqrt[3]*(-(3/2) + 2*t)), 1/2 <= t <= 1}}, 0], 
    Piecewise[{{2^(3*Sqrt[5]*(7/6 - 2*t + (-1 + 2*t)^2)), 
       1/2 <= t <= 1}}, 0]}];




left = Evaluate[Integrate[F[t], t]];

CC = Array[c, {6, 6}];

right = CC . F[t]

I think we will solve left==right for some values of $t$.

$\endgroup$
1
  • 1
    $\begingroup$ How do you measure error? FindMinimum[error[left, right], Flatten@CC] might do it. $\endgroup$
    – Michael E2
    Oct 3 at 15:23
1
$\begingroup$

Using the $L^2$ norm, the error becomes the integral over $t$ of a polynomial in the $c_{ij}$ whose coefficients are functions of $t$. Thus the $L^2$ norm of the error may be computed by numerically integrating the coefficients of the error. This was faster than integrating the error directly or numerically integrating it at each step within the objective function.

vars = Flatten@Array[c, {6, 6}];
errca = MapIndexed[Map[PiecewiseExpand, ##] &, 
   CoefficientArrays[((left - right)^2), vars] (* left, right as in OP *)
   ];
errcai = ArrayRules /@ errca /. (
     # -> 
        Quiet[NIntegrate[#, {t, 0, 1/2, 1}], {NIntegrate::izero}] & /@
       DeleteDuplicates@Flatten@Through[errca@"NonzeroValues"]
     ) // Map[SparseArray];
err = Total@Fold[#2 + # . vars &, Reverse@errcai];
cc = Partition[
   Chop[
    FindArgMin[err, vars],
    1*^-8],
   6];
cc // MatrixForm
$\endgroup$
1
  • $\begingroup$ Thank you. But I couldn't fully understand. Could you explain the mathematical framework of the steps in your code? $\endgroup$
    – 1_student
    Oct 4 at 8:11
0
$\begingroup$

All integrations can also be done analytically in about 70 seconds. This makes it easier to see what is done at the single steps of the whole calculation.

F[t_] = {Piecewise[{{Sqrt[2], 0 <= t <= 1/2}}, 0], 
Piecewise[{{2^(Sqrt[3]*(-(1/2) + 2*t)), 0 <= t <= 1/2}}, 0], 
Piecewise[{{2^(3*Sqrt[5]*(1/6 - 2*t + 4*t^2)), 0 <= t <= 1/2}}, 0],
Piecewise[{{Sqrt[2], 1/2 <= t <= 1}}, 0], 
Piecewise[{{2^(Sqrt[3]*(-(3/2) + 2*t)), 1/2 <= t <= 1}}, 0], 
Piecewise[{{2^(3*Sqrt[5]*(7/6 - 2*t + (-1 + 2*t)^2)), 
  1/2 <= t <= 1}}, 0]};

Plot[Evaluate[F[t]], {t, 0, 1}, 
  PlotStyle -> {Hue[0], Hue[.2], Hue[.4], Hue[5.5], Hue[.7], Hue[.85]}]

(int[t_] = 
Integrate[F[ts], {ts, 0, t}, Assumptions -> 0 < t < 1] // 
 Simplify[#, 0 < t < 1] &) // TableForm // Timing

vars = Flatten[
  Table[ToExpression[StringJoin["c", ToString[i], ToString[j]]], {i, 
6}, {j, 6}], 1]

varsdef = 
  Flatten[Table[
 ToExpression[
StringJoin["c", ToString[i], ToString[j], 
 ToString[_?NumericQ]]], {i, 6}, {j, 6}], 1]

right[t_] = Partition[vars, 6].F[t]

(ii2 = Integrate[
 PiecewiseExpand[(int[t] - right[t])^2 // Expand], {t, 0, 
  1}]); // Timing

(*   {69.891, Null}   *)

tii2[Sequence @@ varsdef] = Total@(ii2);

(famin = FindArgMin[tii2[Sequence @@ vars], vars]) // Timing

(*   {2.796, {-0.118633, 0.567828, -0.0762009, 0.5, 
2.72374*10^-8, -2.69813*10^-8, -0.161573, 0.41647, 2.57896*10^-8, 
0.375172, 9.91451*10^-9, -1.13425*10^-8, -0.042991, 
0.35647, -0.0478373, 0.376858, 
1.47821*10^-8, -3.43372*10^-9, -8.81958*10^-9, -3.11096*10^-9, 
6.75703*10^-9, -0.118633, 0.567828, -0.0762008, -1.54488*10^-8, 
7.88121*10^-9, 3.78811*10^-9, -0.161573, 0.41647, 
1.27111*10^-8, -1.83473*10^-8, 8.6227*10^-9, 
6.70863*10^-9, -0.042991, 0.35647, -0.0478373}}   *)

This is the same result as @Micheal E2 got.

matME2 = {
{-.118633, .567828, -.0762008, .5, 0, 0},
{-.161573, .41647, 0, .375172, 0, 0},
{-.042991, .35647, -.0478373, .376858, 0, 0},
{0, 0, 0, -.118633, .567828, -.0762008},
{0, 0, 0, -.161573, .41647, 0},
{0, 0, 0, -.042991, .35647, -.0478373}
};
solME2 = Flatten[Thread /@ Thread[Partition[vars, 6] -> matME2], 1];

Plot[Evaluate[int[t] - right[t] /. Thread[vars -> famin]], {t, 0, 1}, 
PlotStyle -> {Hue[0], Hue[.2], Hue[.4], Hue[5.5], Hue[.7], Hue[.85]},
PlotRange -> All]

Plot[Evaluate[int[t] - right[t] /. solME2], {t, 0, 1}, 
 PlotStyle -> {Hue[0], Hue[.2], Hue[.4], Hue[5.5], Hue[.7], Hue[.85]},
PlotRange -> All]

Plot[Evaluate[right[t] /. Thread[vars -> famin]], {t, 0, 1}, 
PlotStyle -> {Hue[0], Hue[.2], Hue[.4], Hue[5.5], Hue[.7], Hue[.85]},
PlotRange -> All]

Plot[Evaluate[int[t]], {t, 0, 1}, 
PlotStyle -> {Hue[0], Hue[.2], Hue[.4], Hue[5.5], Hue[.7], Hue[.85]},
PlotRange -> All]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.