Model involving SOD of two variables

  • Context: Graduate 
  • Thread starter Thread starter lol_nl
  • Start date Start date
  • Tags Tags
    Model Variables
Click For Summary

Discussion Overview

The discussion revolves around modeling the velocity of a river stream as a function of its cross-sectional coordinates (y, z). Participants explore the mathematical formulation of the problem, including the derivation of differential equations and boundary conditions, and seek methods to integrate these into a comprehensive function v(y, z) that describes the velocity throughout the river's cross-section.

Discussion Character

  • Technical explanation
  • Mathematical reasoning
  • Homework-related
  • Debate/contested

Main Points Raised

  • One participant describes the modeling of river velocity v(y, z) based on measurements of depth and velocity at various points, leading to a set of differential equations.
  • The participant presents boundary conditions for the problem, including conditions at the riverbanks and the bottom, and discusses the integration of the resulting equations.
  • Another participant points out a potential misinterpretation regarding the constants C[1] and C[2] in the Mathematica code, clarifying that they are function names rather than constants.
  • A later reply acknowledges the clarification and expresses uncertainty about how to solve for C[1] and C[2] using boundary conditions, comparing the problem to the wave equation.
  • Another participant suggests a numerical approach to solving the problem, proposing a simplified version of the differential equation and discussing the setup of a grid for numerical solutions.
  • They share initial results from their numerical method, noting that the derivative with respect to z is close to the expected boundary condition but not exactly zero.

Areas of Agreement / Disagreement

Participants express various viewpoints on the interpretation of the equations and the methods for finding solutions. There is no consensus on the best approach to solve for the functions C[1] and C[2] or on the effectiveness of the proposed numerical methods.

Contextual Notes

Participants acknowledge limitations in their understanding of differential equations and the application of boundary conditions, which may affect their ability to derive solutions. The discussion includes attempts to clarify mathematical interpretations and numerical methods without resolving the underlying complexities.

lol_nl
Messages
40
Reaction score
0

Homework Statement


I'm going to model the velocity of a river stream v(y,z) as a function of the location (y,z) in the cross-section of the river. Here y represents the distance from one of the shores and z the depth of the river.
I have done measurements of the depth z on different locations y, so I can find an approximation for a function z(y) and hence find the profile of the river.
I have also done measurements of the velocity on different points in the river, so now I wish to find a function v(y,z) that describes the speed at every point in the cross-section. This will give me a way to find the average speed by integrating over the function for the cross-area.

Homework Equations


Basically, through a number of assumptions and approximations, the problem can be reduced to a set of two differential equations:

\eta_{y}\frac{\partial^{2} v(y,z)}{\partial y^2} + \eta_{z}\frac{\partial^{2} v(y,z)}{\partial z^2} = -\rho g Sin(\alpha)

By assuming by turns \eta_{y} = 0 and \eta_{z} = 0 I can find functions v(y) and v(z) that describe a kind of an average velocity for any given y or given z.
The problem now is the integration of these two into one function v(y,z).

Given boundary conditions:
* v(0,z) = 0
* v(W,z) = 0 -> W represents max width
* v(y,0) = 0 -> (y,0) is a point on the bottom
* \frac{\partial v}{\partial z}[H] = 0 -> change of speed is zero in for the maximum height H.

The Attempt at a Solution


I've found solved the equations for the assumptions \eta_{y} = 0 and \eta_{z} = 0 and using the boundary conditions found:

v(y) = \frac{\rho g Sin[\alpha]}{2\eta_{y}} (W-y)y

v(z) = \frac{\rho g Sin[\alpha]}{2\eta_{z}} (2H-z)z

The question now is whether there is an easier way to find v(y,z) by directly solving it mathematically, possibly by combining solutions for v(y) and v(z)? I can see that a linear combination of the two would work, but I can't see how I can find the general solution of the function.

I have tried to find a solution using Mathematica, but by filling in the boundary conditions it reduced to the function v(y) that is independent of z. Here is what my code in Mathematica more or less looks like:

Code:
DSolve[{a*D[f[y, z], {y, 2}] + b *D[f[y, z], {z, 2}] == -\[Rho]*g*
    Sin[\[Alpha]]}, f, {y, z}]

{{f -> Function[{y, z}, -((g y^2 \[Rho] Sin[\[Alpha]])/(2 a)) + 
     C[1][(Sqrt[-a b] y)/a + z] + C[2][-((Sqrt[-a b] y)/a) + z]]}}

v[y_, z_] := -((g y^2 \[Rho] Sin[\[Alpha]])/(2 a)) + 
   c1 ((Sqrt[-a b] y)/a + z) + c2 (-((Sqrt[-a b] y)/a) + z);
r1 := v[0, z];
r2 := v[W, z];
r3 := v[y, 0];
r4 := D[v[y, z], z];
{r1, r2, r3, r4}

{c1 z + c2 z, 
 c2 (-((Sqrt[-a b] W)/a) + z) + c1 ((Sqrt[-a b] W)/a + z) - (
  g W^2 \[Rho] Sin[\[Alpha]])/(2 a), (Sqrt[-a b] c1 y)/a - (Sqrt[-a b] c2 y)/
  a - (g y^2 \[Rho] Sin[\[Alpha]])/(2 a), c1 + c2}

Clear[c1, c2]
Solve[{r1 == 0, r2 == 0}, {c1, c2}] 
Solve[{r3 == 0, r4 == 0}, {c1, c2}]
Solve[{r1 == 0, r3 == 0}, {c1, c2}]
Solve[{r2 == 0, r3 == 0}, {c1, c2}]

{{c1 -> (g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]), 
  c2 -> -((g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]))}}

{{c1 -> (g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]), 
  c2 -> -((g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]))}}

{{c1 -> (g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]), 
  c2 -> -((g y \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]))}}

{{c1 -> -((-b g W^2 \[Rho] Sin[\[Alpha]] + b g W y \[Rho] Sin[\[Alpha]] + 
     Sqrt[-a b] g y z \[Rho] Sin[\[Alpha]])/(4 a b z)), 
  c2 -> -((-Sqrt[-a b] g W^2 \[Rho] Sin[\[Alpha]] + 
     Sqrt[-a b] g W y \[Rho] Sin[\[Alpha]] + a g y z \[Rho] Sin[\[Alpha]])/(
    4 a Sqrt[-a b] z))}}

Clear[c1, c2]
c1 := (g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]);
c2 := -((g W \[Rho] Sin[\[Alpha]])/(4 Sqrt[-a b]));
Simplify[v[y, z]]

(g (W - y) y \[Rho] Sin[\[Alpha]])/(2 a)

c1 := -((-b g W^2 \[Rho] Sin[\[Alpha]] + b g W y \[Rho] Sin[\[Alpha]] + 
    Sqrt[-a b] g y z \[Rho] Sin[\[Alpha]])/(4 a b z));
c2 := -((-Sqrt[-a b] g W^2 \[Rho] Sin[\[Alpha]] + 
    Sqrt[-a b] g W y \[Rho] Sin[\[Alpha]] + a g y z \[Rho] Sin[\[Alpha]])/(
   4 a Sqrt[-a b] z));
Simplify[v[y, z]]

(g W (W - y) \[Rho] Sin[\[Alpha]])/(2 a)

In both cases I obtain exactly the solution for v(y).
 
Physics news on Phys.org
Have not gone through it in detail but the first thing I notice is that you seem to be mis-interpreting the line:

\left\{\left\{f\to \text{Function}\left[\{y,z\},-\frac{g y^2 \rho \text{Sin}[\alpha ]}{2 a}+C[1]\left[\frac{\sqrt{-a b} y}{a}+z\right]+C[2]\left[-\frac{\sqrt{-a b} y}{a}+z\right]\right]\right\}\right\}

The C[1] and C[2] are function names and not constants c1 and c2 multiplied by those values in brackets. You realize that or no? So the first one is an arbritrary function of a single variable named C1 at the value given in brackets. Or am I insulting your intelligence by saying that? Sorry if so.
 
Ah, thanks, I actually didn't know that! Indeed, now the solution seems to make more sense to me. It's identical to the wave function except that it is non-homogeneous. Now I presume I can find the functions C[1] and C[2] in the same way as the wave function and then add the factor -\frac{gy^2\rho Sin[\alpha]}{2\eta_{y}} to it?

EDIT: I have simply no idea how to solve for the functions C[1] and C[2] using the boundary conditions. I looked up the solution to the wave equation but the approach there is quite different since it is possible to use the boundary conditions of amplitude and speed at t=0 for x=0. I lack knowledge of differential equations to solve for the functions. Also, if I try to insert the boundary conditions in DSolve in Mathematica, it doesn't give me any solution.
 
Last edited:
lol_nl said:
EDIT: I have simply no idea how to solve for the functions C[1] and C[2] using the boundary conditions. I looked up the solution to the wave equation but the approach there is quite different since it is possible to use the boundary conditions of amplitude and speed at t=0 for x=0. I lack knowledge of differential equations to solve for the functions. Also, if I try to insert the boundary conditions in DSolve in Mathematica, it doesn't give me any solution.

Yeah, me neither. How about solving it numerically? First consider a simplified version:

\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}=-1/2,\quad 0\leq y\leq W,\quad 0\leq z\leq H

\begin{array}{cc} v(0,z)=0 &amp; v(W,z)=0 \\<br /> v(y,0)=0 &amp; \frac{\partial v}{\partial z}\biggr|_{z=H}=0<br /> \end{array}<br />

Or something simple like that. So I'll set up a grid say for now 10x10, do all the difference equations on it with a spacing of one, then solve for all the internal points so that's what 64 linear equations in 64 variables. Alright, I did this quick so not entirely confident of the code but the back-substitution on the results using a sample point in the grid are close to -1/2 which is the first numeric result and the second result is the derivative with respect to z along the top edge at a select point. It's -0.07. Not zero but pretty close to what we want. So then maybe the code is workable so if you want, you can start adding complexity to the PDE and make the grid spacing smaller. When you're done with that, it's not then hard to integrate the interpolation function as you indicated above.

Code:
In[837]:=
Remove[v]
maxv = 10; 

myary = Array[Subscript[v, #1, #2] & , {maxv, maxv}]; 

For[j = 1, j <= maxv, j++, {Subscript[v, 1, j] = 0; Subscript[v, j, 1] = 0; 
     Subscript[v, j, maxv] = 0; Subscript[v, maxv, j] = Subscript[v, maxv - 2, j]; }]; 

myeqns = Flatten[Table[Subscript[v, n, j + 1] + Subscript[v, n, j - 1] + 
       Subscript[v, n - 1, j] + Subscript[v, n + 1, j] - 4*Subscript[v, n, j] == -0.5, 
     {n, 2, maxv - 1}, {j, 2, maxv - 1}]]; 

myvars = Flatten[Table[Subscript[v, n, j], {n, 2, maxv - 1}, {j, 2, maxv - 1}]]; 

mysol = First[Solve[myeqns, myvars]]; 

For[n = 2, n <= maxv - 1, n++, For[j = 2, j <= maxv - 1, j++, 
     Subscript[v, n, j] = Subscript[v, n, j] /. mysol; ]; ]; 

mydatapoints = Flatten[Table[{j, n, myary[[n,j]]}, {n, 1, maxv}, {j, 1, maxv}], 1]; 

they = 2.6; 
thez = 3.8; 
thev = Interpolation[mydatapoints]; 
vy2[y_, z_] = D[thev[y, z], {y, 2}]
vz2[y_, z_] = D[thev[y, z], {z, 2}]
N[vy2[they, thez] + vz2[they, thez]]

vz1[y_, z_] = D[thev[y, z], z]
vz1[6.5, 10]

Out[849]=
InterpolatingFunction[][y, z]

Out[850]=
InterpolatingFunction[][y, z]

Out[851]=
-0.5060005012005363

Out[852]=
InterpolatingFunction[][y, z]

Out[853]=
-0.07272665510909895

I can do a plot of the solution via:

Plot3D[thev[x,y],{x,1,10},{y,1,10}]

which is below. Note at the top edge with respect to z, it's almost flat which I assume is the first partial set to zero there.
 

Attachments

  • myriverplot.jpg
    myriverplot.jpg
    21.8 KB · Views: 470
Last edited:
Thanks a lot. The graph indeed shows the expected result, with the velocity highest midway between the shores and decreasing in the neighbourhood of the shores and bottom.
I'm going to complete it by adding extending the formula and making the grid larger. Also, since I have real data for a number of points (around 20 I believe), I could insert those in my grid and then try to find an Interpolation through all the points. I'll show some results if I have some more.
 

Similar threads

  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 5 ·
Replies
5
Views
3K
  • · Replies 13 ·
Replies
13
Views
3K
  • · Replies 4 ·
Replies
4
Views
3K
  • · Replies 2 ·
Replies
2
Views
4K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 0 ·
Replies
0
Views
3K
  • · Replies 2 ·
Replies
2
Views
3K
  • · Replies 34 ·
2
Replies
34
Views
3K