Anonymous
Anonymous asked in 電腦與網際網路程式設計 · 10 years ago

# matlab程式 求救

y1=x；

y2=x^3-x

Update:

Rating
• Elisha
Lv 6
10 years ago

1. Bisection(二分法)

function r = bisect(fun,xb,verbose)

xeps = 1e-6; % Smallest tolerances are 5*eps

feps = 1e-6;

a = xb(1,1);

b = xb(1,2); % Use first row if xb is a matrix

xref = abs(b - a); % Use initial bracket in convergence test

fa = feval(fun,a);

fb = feval(fun,b);

fref = max([abs(fa) abs(fb)]); % Use max f in convergence test

if sign(fa)==sign(fb) % Verify sign change in the interval

error(sprintf('Root not bracketed by [%f, %f]',a,b));

end

k = 0; maxit = 50; % Current and max number of iterations

while k < maxit

k = k + 1;

dx = b - a;

xm = a + 0.5*dx; % Minimize roundoff in computing the midpoint

fm = feval(fun,xm);

if verbose, fprintf('%4d %12.4e %12.4e\n',k,xm,fm); end

if (abs(fm)/fref < feps) | (abs(dx)/xref < xeps) % True when root is found

r = xm; return;

end

if sign(fm)==sign(fa)

a = xm; fa = fm; % Root lies in interval [xm,b], replace a and fa

else

b = xm; fb = fm; % Root lies in interval [a,xm], replace b and fb

end

end

warning(sprintf('root not within tolerance after %d iterations\n',k));

2. secant(割線法)

function x = secant(fun,a,b)

% where a, b are initial guess

xk = a;

xkm = b;

fk = feval(fun, a);

fkm = feval(fun, b);

kmaxit = 10;

for n = 1:kmaxit

x = xk - fk*(xk - xkm)/(fk - fkm);

f = feval(fun, x);

if (fk - fkm) < 1e-15, break; end

xkm = xk;

xk = x;

fkm = fk;

fk = f;

end

3. function

function f = func(x)

f(1) = x;

f(2) = x.^3-x;

4. main test

% main test

clear all

clc

x1 = bisect(@func,[1,1.5],0)

x2 = secant(@func,1,1.5);

-------------------------------------------------------------------------

我求出的結果,

bisect == > x = 1.0000

secant == > x = 0.9336

您參考看看