Varying Number of Constraints in a Solve Block (Snell's Law, Mathcad 15)

[MoeSzyslak]

Hi, All --

I'm trying to build a solve block with the following constraints:

s1*sin(q1) = s2*sin(q2) = ... = sn*sin(qn) 

However, the number of terms (n) is not know a priori. Is there a way to do this in MC15?

Thanks.

Matt

[LucMeekes]

Is it more difficult than:

Success!
Luc

[MoeSzyslak]

Thanks, Luc -- but that doesn't seem to work -- unless I am doing something else wrong, of course. Please see attached.

[Fred_Kohlhepp]

---

wrote:

Hi, All --

I'm trying to build a solve block with the following constraints:

s1*sin(q1) = s2*sin(q2) = ... = sn*sin(qn) 

However, the number of terms (n) is not know a priori. Is there a way to do this in MC15?

Thanks.

Matt

---

You haven't supplied enough information.  If you're after the s values (s1, s2, etc) and you know the q's, then Luc has given you an answer.  If you know the s's and are looking for the q's, you can use his answer but's going to take some effort.  If you're after both the s's and q's you need more constraints. . .

[MoeSzyslak]

Fred -- you're absolutely correct -- I just posted a sheet in response to Luc's solution. I would appreciate it if you could have a look.

Matt

[Fred_Kohlhepp]

Hmmm

[MoeSzyslak]

EDIT: Added picture of problem geometry.

Thanks, Fred -- I should have been more specific. This is an iterative calculation where I am attempting to find a vector of angles (theta) that is compatible with a vector of slownesses (S). I originally tried the approach you are showing, but it does not converge (see attached). Note that, in reality, I would wrap this whole thing up in a loop with a convergence trigger to jump out, but for just looking at how the first few iterations of the calculations progress, I am just copy-pasting.

Thanks again.

[Werner_E]

Haven't looked at the last sheet you posted and just tried to make the solve block in your first sheet working.

Not sure if the following would do the job. I am not proud of this crude solution and hope someone is able to find a more elegant solution.

[MoeSzyslak]

Thanks for having a look Werner -- pretty slick way to get the solve block working. Unfortunately, the solution you obtain is quite far from the true solution, shown below. What's odd is that the first set of constraints (the sum) is well-satisfied, but the second (the equality) is not.

[MoeSzyslak]

OK, in looking at this again, I see the issue. The ray path lengths are a function of the angles, which will in turn affect the calculation of the slownesses, but this isn't reflected in the solve block. I tried playing around a bit to make it work, but to no avail....

[LucMeekes]

So is this essentially your problem description (where I'm unsure about what are the knowns):

?

Luc

[MoeSzyslak]

Yes, that is the problem. The slownesses (velocities) and the refraction angles are unknown. The vertical distance between measurements, the vector of travel times, and the horizontal offset at the surface are all known.

[Werner_E]

As far as I can see all constraints are pretty much fulfilled to a certain degree.

Its surprising how quick the solve block solves and why changing the values of TOL and CTOL don't noticeably affect the result. Furthermore is Levenberg-Marquardt the only algorithm which finds a solution, the other two (conjugate gradients and Quasi-Newton) fail.

Which are the additional constraints you'd like to add?

[MoeSzyslak]

Yes, that's true, but the matrix of travel lengths (L) is a function of the computed refraction angles. So, we start with an initial guess of slownesses (S=L^-1*T) and use those to compute the refraction angles. The new refraction angles are then used to compute a new travel distance matrix, say L_1. This is used to compute a new slowness vector, say S_1. The calculation is stopped when (S_i+1 - S_i)/S_i < eps where eps is some user-defined convergence criterion.

I didn't include all of this (^^^) in my original solve block, foolishly. Do you think it would be possible to end the solve block with a function and then wrap that function up in a convergence-checking loop?

[Werner_E]

---

Do you think it would be possible to end the solve block with a function and then wrap that function up in a convergence-checking loop?

---

I am not sure cause it was quite unsatisfactory how I rebuilt your third constraint (stacking up T 45 times to have something to compare with.

After all when using a solve block you should not think in iterations and how one matrix is constructed using the result of the aforegoing iteration step. You should think in results, in equations, in constraints you want to be fulfilled. You want to turn the solve block into a function (of what) and use it in a self written iteration? Why won't you let do the solve block the whole job?

[MoeSzyslak]

Unfortunately, I don't think it is possible to let the solve block do the whole job -- please have a look at the last worksheet I posted, which outlines the calculation procedure step-by-step.

[Werner_E]

---

wrote:

Unfortunately, I don't think it is possible to let the solve block do the whole job -- please have a look at the last worksheet I posted, which outlines the calculation procedure step-by-step.

---

I must confess that I don't see much similarity in that new sheet compared to the one you posted first.

All I see are two solve blocks and only value V2 from the first solve block is used in the next one and V2 is not changed. The angles you calculated in the first solve block are not used again, correct?

Is it a typo that you have V2 (and not V3) in the denominator of the last term in your first constraint?

I also don't see how and why the third (not working) solve block would do the very same as the first two!?

[MoeSzyslak]

OK, I was trying to take some shortcuts with this, but perhaps it's best to fully explain the problem. Attached I have produced the iterative solution (brute force) for the first three layers. The solution is, as best I can tell, correct. However, I am looking for a way to formulate it such that I can pass it vectors of layer thicknesses and travel times of arbitrary length -- clear the attached approach is not amenable to that.

Thanks again, all.

[Fred_Kohlhepp]

But you're changing the problem each time.

• You solve for the velocity of the first layer:  V1 = sqrt(Z1^2+S^2)/T1  Implying that Z1 is the horizontal distance for a 3 meter thick layer.
• Then your first solve block says that the sum of z tan(theta) equals the vertical height.  But all of your Z's are 1 meter  (Should be z* cot(theta)   ?)  So now you have the same three meter vertical thickness and two meters horizontal travel.  Not the same problem!

[MoeSzyslak]

Hi, Fred et al. --

The layer thickness is 1 m, but the horizontal offset of the signal source is 3 m. In general, of course, these values are arbitrary.

I have attached an image of the problem geometry. The problem isn't changing, I am just not going a good job explaining it. I had hoped that I could solve it in one way when I started this thread, but now I am pretty convinced that way won't work. So I have re-cast it a couple of different ways in hopes of better explaining it, but I have served only to sow confusion.

In the attached figure, We have a wave source on the surface and some receivers on a vertical line some horizontal distance away, S. The travel times from source to receivers are measured. In the first layer (V_1), the wave velocity (or slowness) may be calculated directly. That velocity is then used to calculate the wave velocity in layer two (V_2), which requires the knowledge of the refraction angles q_12 and q_22. Thus, the problem is now iterative. Once V_2 is known, the calculation is continued for subsequent layers. In the figure, X_ij corresponds to a property of the jth raypath in the ith layer where X is some quantity (length, angle, etc.).

Does this make sense?

Thanks.

Matt

[Werner_E]

> Does this make sense?

Yes, it does. Now I think I know what you are after and understand the equations and the approach you had used in your second sheet.

Unfortunately this does not mean that I would have a solution at the moment.

[Fred_Kohlhepp]

I understand now.

The good news is that I can "automate" it.

The bad news is that it doesn't give the right answers.

[Werner_E · Accepted solution]

OK, I couldn't resits 😉

Here's my approach. Looks like the values are close enough to the one you expected and all conditions should be fulfilled. You may even use "Find" instead of "MinErr".

The use of the predefined guess vector for the angles isn't really elegant but the alternative would be to create a suitable guess vector of the correct size before calling the solve block function and add it as additional function argument.

I guess you are only interested in the velocities and not in the various values of the angles in each iteration step. If you should need all the angles, it should be easy enough to collect them in a matrix and return them as well.

BTW, is it possible that the time values you posted are a bit inaccurate 😉

[MoeSzyslak]

This is fantastic, Werner -- very ingenious.

Yes, the time values were imperfect -- just some synthetic numbers to test the algorithm. I'm going to give it a go with some real data.

On a related note, there do exist algorithms that accept not just absolute travel times, but also standard deviations. These are true inversions with regularization, obviously, Those approaches are described in the literature often in a "hand wavy" fashion, making them difficult to implement if you don't want to start at the beginning....