MATLAB Answers

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
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])
%
% See also: fsolve (in optimization toolbox)
% Author: Are Mjaavatten, Telemark University College, Norway
%
% The update formula for the Jacobian J is taken from
% the broyden function written by Matthias Heinkenschloss:
% http://read.pudn.com/downloads158/sourcecode/math/705810/sysnlineq/broyden.m__.htm
%
% 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
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
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
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?

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

回答 (0 件)

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by