# This is Broyden.m I found. I have a system of 3 functions and 6 unknowns, I think this function is meant to solve 1 unknown per equation or am I not coding my eqns() function corrrectly?(in comments)

11 ビュー (過去 30 日間)
Timothy Nsubuga 2019 年 7 月 1 日
コメント済み: Timothy Nsubuga 2019 年 7 月 1 日
function [x, ithist] = broyden(f,x0,opt,bounds)
% Solve system of nonlinear equations by Broyden's method
% Broyden's method is a quasi Newton method which updates an approximate
% Jacobian at each new Newton step.
%
% Usage:
% [x, ithist] = broyden( f, x0, opt, bounds)
%
% Input parameters:
% f name of a matlab function that evaluates F. (Required)
% x0 iteration starting point (Required)
% opt Options struct. Fields: (Optional)
% tolfun Function tolerance. (Default tolfun = 1.e-10)
% Iterations end if ||F(x)||_2 < tolfun
% tolx Step length tolerance (Default tolx = 1.e-8)
% Newton's method stops if ||step||_2 < tolx,
% where step is the Newton step.
% maxiter maximum number of iterations (Default maxiter = 50)
% bounds Optional bounds on x Default: no bounds.
% bounds(:,1) <= x < bounds(:,2)
%
% Output parameters:
% x approximation of the solution.
% ithist iteration history struct. Fields: x, f(x), norm(f)
%
% Example:
% f = @(x) [x(1)^2 + x(2)^2 - 4; exp(x(1)) + x(2) - 1];
% x = broyden(f,[1;1])
% The other root has x(1) > 0, x(2) < 0. This may be found by selecting
% another inital x, or by using constraints. E.g.:
% x = broyden(f,[1;1],[],[0 2.5;-2.5 0])
%
% Author: Are Mjaavatten, Telemark University College, Norway
%
% The update formula for the Jacobian J is taken from
% the broyden function written by Matthias Heinkenschloss:
%
% Version 1: 2015-12-29
% Version 1.1: 2016-11-29 Added iteration history
optfields = {'maxiter','tolfun','tolx'};
defaults = {50,1e-10,1e-8,[]};
if nargin < 3
opt = [];
end
for i = 1:3
if isfield(opt,optfields{i})
if isempty(opt.(optfields{i}))
opt.(optfields{i}) = defaults{i};
end
else
opt.(optfields{i}) = defaults{i};
end
end
x = x0(:);
it = 0;
F = feval(f, x);
if ~(size(x,1) == size(F,1))
error('f must return a column vector of the same size as x0')
end
normf = norm(F);
J = jacobi(f,F,x); % Intial Jacobian matrix
if nargout > 1
ithist.x = [x(:)';zeros(opt.maxiter,length(x))];
ithist.f = [F(:)';zeros(opt.maxiter,length(x))];
ithist.normf = [normf;zeros(opt.maxiter,1)];
end
normdx = 2*opt.tolx;
while( it < opt.maxiter+1 && normdx > opt.tolx && normf > opt.tolfun)
if rcond(J) < 1e-15
error('Singular jacobian at iteration %d\n',it)
end
dx = -J\F;
normdx = norm(dx);
if nargin > 3 % variable bounds are supplied
% make sure x stays within bounds
for j = 1:20
jl = find(x+dx<bounds(:,1));
dx(jl) = dx(jl)/2;
ju = find(x+dx>bounds(:,2));
dx(ju) = dx(ju)/2;
if isempty(jl) && isempty(ju)
break
end
end
end
x = x+dx;
it = it+1;
F = feval(f, x);
normf = norm(F);
J = J + F*dx'/(dx'*dx);
if nargout > 1
ithist.x(it+1,:) = x(:)';
ithist.f(it+1,:) = F(:)';
ithist.normf(it+1,:) = normf;
end
end
% Check if the iterations converged and issue warning if needed
if it >= opt.maxiter && norm(F) > opt.tolfun
warning('No convergence in %d iterations.\n',it+1)
elseif normf>opt.tolfun
warning('Newton step < %g, but function norm > %g\n',...
opt.tolx,opt.tolfun)
elseif normdx>opt.tolx
warning('Function norm < %g, but newton step norm > %g\n',...
opt.tolfun,opt.tolx)
end
if nargout > 1
ithist.x(it+2:end,:) = [];
ithist.f(it+2:end,:) = [];
ithist.normf(it+2:end) = [];
end
end
function J = jacobi(f,y0,x)
% Quick and dirty numerical Jacobian for function f at x
% y0: f(x);
delta = 1e-6*(max(1,sqrt(norm(x))));
n = length(y0);
m = length(x);
J = zeros(n,m);
for i = 1:m
dx = zeros(m,1);
dx(i) = delta/2;
J(:,i) = (feval(f,x+dx)-feval(f,x-dx))/delta;
end
end

#### 3 件のコメント

Timothy Nsubuga 2019 年 7 月 1 日
%*************************************************************************
%%This fucntion will store transducer coordinates and differences of
%%arrival times between each transducer and transducer 1,
%the arrival times are aquired via graphData() and is an input broyden() to solve for the
% signal location and travel time to transducer 1, adding the difference in
% arrival time of each transducer to this travel time yields the travel
% time
%%
%**************************************************************************
function [fun] = eqns(u)%unknown values are input parameters
%the pop rock falls vertically and so xs = zs = 0 and we assume the center
%of the borehole is at y = 0
xs = 0;
zs = 0;
%Transducer 1 and 2 locations in meters based on 3D print sample with a
%fracture 1/3 of the sample height
xt1 = -0.05;
yt1 = 0.015;
zt1 = -0.015;
xt2 = -0.05;
yt2 = 0.015;
zt2 = 0.015;
xt3 = -0.05;
yt3 = -.005;
zt3 = -0.015;
v = 2500;% wave speed in material determined by AST and distance formula
%differnce in arrival time between tranducer 1 and transducer 2;
% acquire locations of first peaks in transducers 1 and 2
File_names = dir('*.csv'); %attain every waveform in the event folder
%get files of interest, in this case those of transducer 1 and 2
file1 = File_names(3).name;
file2 = File_names(4).name;
file2 = File_names(5).name;
plus1 = 1; %offset y values
plus2 = 2;
plus3 = 3;
hold on
% %prof requested the dts so I need all the first peak locations
% numFiles = length(File_names);
% incr = 1:numFiles;
% for i = 1:numFiles %call graphData() for each waveform
% file = File_names(i).name;
% hold on
% plus = incr(i);%needed to offset graphs for plotting
% [p,l] = graphData(file,plus);
%
%
% end
[p1,l1] = graphData(file1,plus1);
% (ORIGINAL CODE IS COMMENTED OUT FOR NOW TO GET ALL
% DTS FROM EVENT 1)
% L1 = l1;
[p2,l2] = graphData(file2,plus2);
[p3,l3] = graphData(file2,plus2);
% L2 = l2;
% graphData() gives me locations on the first peak but I only require one of
% them
% maxL1 = max(L1);
%
% maxL2 = max(L2);
%multiply each location by the sample interval, 1*10^-7 to get the arrival
%times for transducers 1 and 2
art1 = l1*1e-7;
art2 = l2*1e-7;
art3 = l3*1e-7;
%
% tran10= 1038*1e-7;
%
% tran13 = 1019*1e-7;
%
% tran1 = 1047*1e-7;
%
% tran2 = 1054*1e-7;
%
% tran3 = 1038*1e-7;
%
% tran5 = 1048*1e-7;
%
% tran6 = 1042*1e-7;
%
% tran8= 1022*1e-7;
%
% tran9 = 1034*1e-7;
%
% dt10 = tran10 - tran1;
% dt13 = tran13 -tran1;
% dt1 = tran1- tran1;
% dt2 = tran2- tran1;
% dt3 = tran3- tran1;
% dt5 = tran5- tran1;
% dt6 = tran6- tran1;
% dt8 = tran8- tran1;
% dt9 = tran9 - tran1;
%deltat1 is clearly zero
deltat2 = art2 - art1;
deltat3 = art3 - art1;
% deltat2 = 7e-7;
ys = u(1); %source location and travel time from source to transducer 1
% are unkown and will be guessed initally with values of
% [0.5;2e-5] as input argumnents to eqns()
ts = u(2);
%non-linear system of equations, the first is for transducer 1 and the
%second is for transducer 2
fun(1,1) = (sqrt((xt1- xs)^2 + (yt1 - ys)^2 + (zt1 - zs)^2)/v )- (ts);
fun(2,1) = (sqrt((xt2- xs)^2 + (yt2 - ys)^2 + (zt2 - zs)^2)/v) - (ts - deltat2);
fun(3,1) = (sqrt((xt3- xs)^2 + (yt3 - ys)^2 + (zt3 - zs)^2)/v) - (ts - deltat3);
end
%end of eqns()*************************************************************
Walter Roberson 2019 年 7 月 1 日
You are right that the code cannot have a different number of variables than it has return values from the function.
It might be interesting to return 3 extra 0 from the function.
Timothy Nsubuga 2019 年 7 月 1 日
guess = [0.5;2e-5];
>> f = @eqns;
>> x = broyden(f,guess);
Here is my call, each equation has ts and ys, and I want to solve for these 2 variables. guess is the corresoponding values I beleive ys and ts are, any tips to get this to work?

サインインしてコメントする。

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by