3
$\begingroup$

I am trying to use the NIntegrate function to find the difference in the area beneath two 2 dimensional splines (how much area of the river bed was lost or gained from one year to the next) but, I cannot figure out how to do this repeatedly at 1m intervals along the splines in a way that will give me values I can add directly to a plot (without having to copy/paste or hard-code the output). An ideal output would resemble something similar to {{1,x1}, {2,x2},etc.} or be a format that can be converted to something plottable (I'm aware the explanation may be confusing but I am teaching myself Mathematica with no coding and barely any useful computer knowledge beyond Microsoft Office.

Below is the only loop method I have been able to use that gives me the value of the NIntegrate of each interval, but the output for each calculation is its own individual cell and obvious Prints it, rather than being in a useful format (as far as I can tell). "Testing" is Spline1 minus Spline2 and the image below is a screenshot of part of the output of this line of code.

For[i = 0, i < 21, i++, Print[NIntegrate[Testing, {x, i, i + 1}]]]

enter image description here

This screenshot is a plot of the 2 splines being compared. enter image description here

Any help or tools I could be directed to that might help solve this issue would be greatly appreciated.

$\endgroup$
1
  • 2
    $\begingroup$ Not sure I understand what output you are looking for. May be this? testing = x^2/3; {#, NIntegrate[testing, {x, #, # + 1}]} & /@ Range[0, 20] This gives {{0, 0.111111}, {1, 0.777778}, {2, 2.11111},....etc. which you can now plot using ListPlot for example $\endgroup$
    – Nasser
    Nov 13, 2022 at 23:03

2 Answers 2

6
$\begingroup$

Given two splines, with simple examples chosen below:

spline1[x_] := 1.21 x^2;
spline2[x_] := (x+1)^2;

Plotting these:

Plot[{spline1[x], spline2[x]}, {x, 0, 10}]

spline1 and spline2

Then the accumulated difference between them over a range (e.g. x from 0 to 10) can be found using NDSolve:

sol = NDSolve[{d'[x] == spline2[x] - spline1[x], d[0] == 0}, d, {x, 0, 10}]

This saves the continuously accumulated difference starting at $x=0$ as a function d[x] in sol. This difference can be plotted by:

Plot[d[x] /. sol, {x, 0, 10}]

accumulated difference between spline2 and spline1

NDSolve is often useful machinery for collecting numerical integration results over a range of input values. Since we know these functions are continuous, it should also be the case that for any x1 and x2, NIntegrate[spline2[x] - spline1[x], {x, x1, x2}] == (d[x2] - d[x1] /. First[sol]).

NIntegrate[spline2[x] - spline1[x], {x, Pi, 2 Pi}]
d[2 Pi] - d[Pi] /. First[sol]

17.5573

17.5573

First is only necessary here because NDSolve can potentially find multiple solutions for much more complicated systems, and d[x] /. sol would return all of them. In this case, there's only the one so the first one is also the one we are interested in.

$\endgroup$
1
  • $\begingroup$ I managed to incorporate parts of the method you used to get my desired result. Thank you for the help. I will try to add my final result below to show exactly what it was I needed. $\endgroup$ Nov 27, 2022 at 3:31
1
$\begingroup$

Below was my adaptation of the answers above to get the result needed (I am the one who initially asked the question)

This is my first time doing any type of coding and teaching myself so some aspects may look slightly primitive and be obvious to some people.

Finding the Difference between the two splines.

DamThree1u2119 = SpDamThree1u2021[x] - SpDamThree1u2019[x]

NIntegrate to get the difference in area beneath each spline at 1m intervals, showing the positive and negative differences between the data taken in 2019 and 2021

DamThree1u2119AreaDif = NIntegrate[DamThree1u2119, {x, #, # + 1}] & /@ Range[0, 30]

giving the below values

 DamThree1u2119AreaDif = {0.549036,0.734469,0.014559,0.0222887,0.0912996,0.0882717,0.100567,0.133802,0.0854523,0.0744483,0.068068,0.119126,0.135407,0.141813,0.113952,0.095985,0.0809483,0.0797658,0.135477,0.056064,-0.0131309,-0.0272682,-0.0428288,-0.0856341,-0.248317,-0.354404,-0.337383,-0.506746,-0.568715,-0.716076,-0.854683}

Code used to get the BarChart shown below

BarChart[DamThree1u2119AreaDif,BarSpacing->0,ColorFunction->Function[{y},If[y>0,Green,Red]],ColorFunctionScaling->False,PlotRange->{-1.5,1.5},Axes->{False,True},ImageSize->300,AxesLabel->{None,HoldForm["Area\nchange\n(m^2)"]},LabelStyle->{GrayLevel[0]},Background->LightGray]

enter image description here

Code used to display the splines to see downcutting/filling (red/green respectively) of a cross profile of a river bed, shown below.

DifDamThree1u2119=Plot[{SpDamThree1u2021[x],SpDamThree1u2019[x]},{x,0,31},PlotRange->{{0,31},{96,101}},Filling->{1->{2}},FillingStyle->{Red,Green},PlotStyle->{Directive[RGBColor[0.77,0.45,0.98],Thick],Directive[RGBColor[0.7,0.69,0.3],Thick]},Background->LightGray,PlotLabel->"1u 2019-2021",AxesLabel->{HoldForm["(m)"],HoldForm[Elevation[m]]},ImageSize->300]

enter image description here

Looking at river bed change between the years 2018, 2019, and 2021 combined

Title1u=Panel[Style["Cross Profile Change 1u",White,20],ImageSize->900,Background->Lighter[Gray, 0.5],Alignment->Center];
Deploy@Grid[{{Title1u,SpanFromLeft},{DifDamThree1u2119,DifDamThree1u1918,DifDamThree1u2118},{DifChartDamThree1u2119,DifChartDamThree1u1918,DifChartDamThree1u2118}},Spacings->{0,0},Dividers->{All,False}]

enter image description here

Thanks to everybody who answered and hope this can be of some use to anybody else.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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