Trying to do an fsolve via an anonymous function

I am trying to compute something from a set of nonlinear algebraic equations. I have a function I called initial which defines the set of nonlinear algebraic equations, I use a function initial.m as my function that I need to set as as f(x)=0. I wanted to try and see how fzero sorts this problem, I would normally code up Newton's method to do this.

The idea behind this is that I'm looking at a rod ofdensity rho_0, and I want to split it up into N equally heavy sections. The problem I have is to find how long those intervals are. I'm approximating the mass as 0.5*(rho_0(X(i))+rho_0(X(i+1)))*(X(i+1)-X(i)). It runs but it gives me zeros and bombs out.

Do I need to do back to coding up Newton's method?

[Steven Lord]

Are you trying to use fsolve from Optimization Toolbox or fzero from MATLAB? The question states fsolve but then you said "I wanted to try and see how fzero sorts this problem" and this is tagged "fzero".

The fzero function states that the function fun you pass in as the first input must "accepts a scalar x and returns a scalar fun(x)." and it sounds like you want to find roots of a system of equations in multiple unknowns.

Also, what does "bombs out" mean in this context?

• Do you receive warning and/or error messages? If so the full and exact text of those messages (all the text displayed in orange and/or red in the Command Window) may be useful in determining what's going on and how to avoid the warning and/or error.
• Does it do something different than what you expected? If so, what did it do and what did you expect it to do?
• Did MATLAB crash? If so please send the crash log file (with a description of what you were running or doing in MATLAB when the crash occured) to Technical Support so we can investigate.

[Mat]

I'm using fsolve from matlab, which I'm told doesn't work with systems. I think I will have to do this manually..

[Stephen23]

"I'm using fsolve from matlab, which I'm told doesn't work with systems."

The FSOLVE documentation states "fun is a function that accepts a vector x and returns a vector F, the nonlinear equations evaluated at x". FSOLVE is intended for solving systems of equations, your information source is wrong.

source: https://www.mathworks.com/help/optim/ug/fsolve.html#buta__s-2

Note that your uploaded code actually uses FZERO (not FSOLVE): mixing up functions like this is unlikely to produce meaningful results. Reading the documentation is more likely to produce meaningful results.

I would not recommend doing this "manually", MATLAB has plenty of optimization tools that work.

[Mat]

I've noticed that although matlab has a great deal of inbuilt functions, it is often easier and better to write your own functions for doing these things. It's often quicker in the long run.

[Steven Lord]

What leads you to believe that "It's often quicker in the long run." to write your own functions instead of using the ones included in MATLAB or other MathWorks products? Are you talking about development and maintenance time, running time, time spent learning how to use the function, or something else? I admit if you write your own function the "time spent learning how to use the function" is likely going to be small -- but how about the time spent relearning it when you try to use the code again in a month or two?

Generally, like in cases such as this one from about two weeks ago, attempting to replicate the full behavior of a built-in function can be tricky. There can be a lot of "devil in the details" cases that you may not think about. Granted, if all you need is a small piece of the functionality implementing that specific behavior may be quicker, but often problems (like gases and cats in boxes) tend to expand to fill more space than you expect.

Answer 1

[Mat]

If you use the trapezium rule to approxumate the integral, then the problem for finding the co-ordinates just becomes a linear problem, and you can just invert a matrix.

[Steven Lord]

you can just invert a matrix.

Please don't do that. If you need to solve a system of linear equations, use the backslash operator instead of inverting and multiplying.

Answer 2

[Star Strider]

I exrtracted your code sections and attempted to run it.

You need to check the trapz integration. The call to it is correct, however you are giving it scalars, not vectors, and it cannot integrate scalars.

My analysis --

% ----- rho_0 -----

function y=rho_0(X,L)
  z=0.5*exp(-(X-0.5*L).^2);
  X
  z
  M=trapz(X,z)
  y=z/M;
 end%function

% ----- initial -----

function y=initial(X,N,dh,L)
%This function
y=zeros(N+1,1);
for i=1:N
    i
    X(i)
    L
    S1 = rho_0(X(i),L)
    S2 = rho_0(X(i+1),L)
  y(i)=0.5*(rho_0(X(i),L)+rho_0(X(i+1),L))-dh
end%for

X(1)=0;
end%function

% ----- new_idea -----

%This code is yet another attempt to solve the 1D sintering equations via a differnt approach.
%The incremental mass dh, will be held constant but the initial co-ordinates X, will be non-uniform.

%The initial co-ordinates are computed as a solution of an systen of non-linear algebraic equations.
clear; clc;

%initial length of rod
L=1;

%Amount of mass
M=1;

%Split the rod into N elements;
N=200;
X_0=linspace(0,L,N+1);

%set dh=constant
h=linspace(0,M,N+1);
dh=h(2);

%initial density as a function of X_0:
rho_start=rho_0(X_0,L);

X = 1×201

0 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 0.0350 0.0400 0.0450 0.0500 0.0550 0.0600 0.0650 0.0700 0.0750 0.0800 0.0850 0.0900 0.0950 0.1000 0.1050 0.1100 0.1150 0.1200 0.1250 0.1300 0.1350 0.1400 0.1450

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

z = 1×201

0.3894 0.3913 0.3933 0.3952 0.3971 0.3990 0.4009 0.4028 0.4046 0.4065 0.4083 0.4102 0.4120 0.4138 0.4156 0.4174 0.4191 0.4209 0.4226 0.4244 0.4261 0.4278 0.4295 0.4311 0.4328 0.4344 0.4360 0.4376 0.4392 0.4408

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

M = 0.4613

%Create an anonymous function from initial
disp("N = ")

N = 

N

N = 200

disp("dh = ")

dh = 

dh

dh = 0.0050

disp("L = ")

L = 

L

L = 1

myfun = @(z) initial(z,N,dh,L);

%Get the initial non-uniform Lagrangian co-ordinates
% X=fzero(myfun,X_0);
X=fsolve(myfun,X_0)

i = 1

ans = 0

L = 1

X = 0

z = 0.3894

Error using matlab.internal.math.getdimarg
Dimension argument must be a positive integer scalar within indexing range.

Error in trapz>getDimArg (line 90)
dim = matlab.internal.math.getdimarg(dim);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in trapz (line 44)
dim = min(ndims(y)+1, getDimArg(dim));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in solution>rho_0 (line 7)
M=trapz(X,z)
^^^^^^^^^^^^
Error in solution>initial (line 20)
S1 = rho_0(X(i),L)
^^^^^^^^^^^^^^^^^^
Error in solution>@(z)initial(z,N,dh,L) (line 61)
myfun = @(z) initial(z,N,dh,L);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fsolve (line 167)
fuser = feval(funfcn{3},x,varargin{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.

.

[Mat]

trapz is fine. I'm simply computing the total mass given the initial density.

[Star Strider]

Calling trapz with scalar arguments will produce no useful information.

Answer 3

[John D'Errico]

You say you want to use fzero. But you need to understand that fzero is a ONE variable solver. It CANNOT solve a multi-variable problem.

Do you need to use Newton's method? No, not at all. In fact, it is rarely ever a good idea to write your own rootfinding codes, certainly not if you are as unsure about the problem as your question seems to imply.

Your question is confusing however. You talk about a muti-variable prob;em, but then you state you are looking to identify points (in ONE dimension) along a rod. And I'd not even use a root finder for that. I'd use interp1, at least if I understand your question. You can visualize the integral of rho_0 as a monotonic function. Use cumtrapz to compute the cumulative integral.

Now interpolate a set of uniformly spaced points in cumulative rho.

[Mat]

I find that coding up Newton's method is the only thing that works in most cases.

I thought fzero works for systems as well as single equations.

[Stephen23]

"I thought fzero works for systems as well as single equations."

The FZERO documentation states that the objective function "fun accepts a scalar x and returns a scalar fun(x)." In general a scalar cannot represent a system of equations without loss of information: FZERO is very unlikely to help with your system of equations.

source: https://www.mathworks.com/help/matlab/ref/fzero.html#btoc6lj-1-fun

Answer 4

[Matt J]

It doesn't look like you need an iterative solver at all. It appears that your rho0() function is just a reinvention of normcdf (perhaps with some unit scaling). The solution will therefore just be some appropriately scaled version of,

X=norminv(  linspace(0,1,N),  L/2)  %closed-form solution

[Mat]

This was simply an example, it would be any continuous function.

[Matt J]

The extension to arbitrary continuous functions is no big deal. You just compute its inverse CDF and proceed similarly.