Корни нелинейного уравнения с использованием Matlab

Я использую Matlab, чтобы найти корни нелинейной функции. Уравнение длинное, и я использовал еще один .m для сохранения функции, код которой выглядит так:

function x_c = f_x_c(s,H,VA,Lo,qc,EA,NF,Sj,Fj)

if (s < 0) || (s > Lo);
    disp('The value of s is invalid')
    disp(['s = ' num2str(s)]);
    return
end

C1 = H/qc;
if NF == 0
    n = 0;
    sn = 0;
    sum_Fj = 0;
end
if NF >= 1
    Sj_Q = [0; Sj; Lo];
    %Determine n and sn if 0 <= s < Lo:
    if s < Lo
    STOP = 0;
    k = 0;
        while STOP == 0
        k = k + 1;
            if (s >= Sj_Q(k,1)) && (s < Sj_Q((k + 1),1))
            STOP = 1;
            end
        end
    n = k - 1;
    sn = Sj_Q(k,1);
    end

    %Determine n and sn if s = Lo:
    if s == Lo
        n = NF;
        sn = Sj(NF,1);
    end

sum_Fj = sum(Fj(1:n,1));
end


x_c = (H/EA)*s;
x_c = x_c + C1*asinh((qc*s - VA + sum_Fj)/H) + ...
- C1*asinh((qc*sn - VA + sum_Fj)/H);

for j = 1:n
    sk = Sj_Q((j + 1),1);
    sk_1 = Sj_Q(j,1);
    sum_Fj = sum(Fj(1:(j - 1)));
    x_c = x_c + ...
    + C1*asinh((qc*sk - VA + sum_Fj)/H) + ...
    - C1*asinh((qc*sk_1 - VA + sum_Fj)/H);
end

Здесь переменная H. С кодом проблем нет, потому что он возвращает мне это длинное уравнение, когда я набираю следующее в основном файле.

syms x
equation = -(XB - XA) + f_x_c(s,x,VA,Lo,qc,EA,NF,Sj,Fj);  %Replaced H with variable H and all other arguments are pre-defined

Теперь я хочу решить это уравнение около H0. Когда я ставлю fzero(@(x)equation, H0), это дает мне ошибку, которая выглядит так

Undefined function 'isfinite' for input arguments of type 'sym'.

Error in fzero (line 308)
    elseif ~isfinite(fx) || ~isreal(fx)

Error in main (line 50)
fzero(@(x)equation, H0)

Как я могу решить эту проблему?

ИЗМЕНИТЬ:

У уравнения есть по крайней мере один корень, потому что, если я использую ezplot для построения графика функции, я получаю следующий рисунок.введите описание изображения  здесь


person India Slaver    schedule 02.12.2014    source источник
comment
@patrik Нет, у меня нет такой переменной. Фактически, в fzero.m есть одна строка, в которой есть isfinite. Эта строка выглядит как if any(~isfinite([fa fb])) || any(~isreal([fa fb]))   -  person India Slaver    schedule 02.12.2014
comment
Извини, я виноват. Однако я попытался воспроизвести ошибку с помощью более простой функции, но в настоящее время не могу воспроизвести ошибку. Строка function x_c = equation = XB - XA + f_x_c(s,x,VA,Lo,qc,EA,NF,Sj,Fj); также не может быть запущена. Могу я спросить вас, что вы хотите сделать? У вас есть файл функции, который принимает входные данные и дает выходные данные, из которых вы хотите найти корень?   -  person patrik    schedule 02.12.2014
comment
@patrik Да, Sj и Fj должны быть матрицами в функции, поэтому ваш код может не работать. Код слишком большой! Вы хотите, чтобы я выложил все это здесь?   -  person India Slaver    schedule 02.12.2014
comment
Нет, скорее нет. Но я имею в виду, что syms x; function x_c = equation = XB - XA + f_x_c(s,x,VA,Lo,qc,EA,NF,Sj,Fj); должно давать ошибку компиляции. В-третьих, действительно ли здесь требуется символьный набор инструментов?   -  person patrik    schedule 02.12.2014
comment
@patrik Если я просто напишу это, код работает нормально. Проблема возникает, когда я использую функцию fzero. Если я использую disp(x_c), это показывает мне то большое уравнение, которое я пытаюсь решить.   -  person India Slaver    schedule 02.12.2014
comment
Но тогда этот syms x; function x_c = equation = XB - XA + f_x_c(s,x,VA,Lo,qc,EA,NF,Sj,Fj); не может быть таким, каким он выглядит, потому что это дает ошибку компиляции. Ключевое слово function так не работает! Кроме того, x_c = equation = XB... не является правильным синтаксисом. Если вам нужна помощь, пожалуйста, дайте пример кода, который работает. У меня есть идея, в чем проблема, но если я точно не знаю, как вы получаете ошибку, я не смогу вам помочь. Я мог бы сказать, что вам нужно переписать код совершенно по-другому, который работает, но я предполагаю, что это не ваш вариант.   -  person patrik    schedule 02.12.2014
comment
@patrik x_c = equation = XB... здесь был ошибкой. Я этого не заметил. Хотя в моем коде это правильно. Есть 5 файлов кода. Мне их всех сюда выложить? Или я могу связаться с вами в другом месте? Пожалуйста, помогите мне!   -  person India Slaver    schedule 02.12.2014
comment
@patrik Пожалуйста, посмотрите на график уравнения. Он определенно пересекает ось x один раз.   -  person India Slaver    schedule 02.12.2014


Ответы (2)


Я вижу, я попал почти в то же место, что и вы, однако я добрался до строки 309, так как, видимо, использование isfinite на sym не вызывает сбоя в matlab 2014b. Однако isfinite возвращает false, по-прежнему выдавая ошибку. Но в точку: кажется, что fzero не должно использоваться для символов. Вы уверены, что символы здесь необходимы? Глядя на вызов функции, кажется, что это не обязательно, поскольку вывод выглядит так, как будто он в любом случае должен быть числовым. Также вы можете утверждать, что

syms x; fzero(@(x) x^2-1,1)

работает, но тогда x в @(x) имеет более высокий приоритет, который не является символом (в программировании мы говорим, что переменная имеет другую область действия).

Если важно, чтобы x было sym, вместо этого следует использовать решатель уравнений solve, который работает для символьного вывода.

syms x;
equation = x^2-3;
solve(equation,'x')

Однако это может привести к сбою или дать вам невероятно длинные ответы для сложных функций (а также для выражений, не имеющих хороших дробных ответов, попробуйте syms x; equation = x^3-3.17+1.37*x;). Это также медленнее и поэтому не рекомендуется, если у вас нет произвольных констант в выражении (например, f = x^2-a ‹=> x = +- a^(1/2), где a должно быть определено позже или вы хотите что-то сделать с решением для нескольких значений a).

person patrik    schedule 03.12.2014

Я не знаю, что вызвало проблему, поскольку я не знаком с символическим пакетом в Matlab.

Однако вы должны быть в состоянии решить это довольно легко, используя какой-нибудь численный решатель. Из того, что я могу прочитать о вашем коде, вы хотите решить приведенное выше уравнение равным нулю, это именно то, что делает Ньютон Рэпсон. Вы можете найти метод здесь: http://en.wikipedia.org/wiki/Newton 's_method

Поскольку вы, вероятно, не знаете точную производную своей функции, вы можете оценить ее, используя приближение первого порядка, где вы просто используете определение дифференциалов, но поскольку мы не можем допустить, чтобы h равнялось 0, мы выбираем h очень маленьким, в Matlab я обычно используя sqrt(eps). Таким образом, приближение становится следующим: f'(x) = (f(x+sqrt(eps))-f(x))/sqrt(eps).

В противном случае вы можете использовать мой метод, который работает в 1 измерении, вы можете найти его здесь: http://imada.sdu.dk/~nmatt11/MM533-2014/ Вам необходимо скачать fpisimple.m и mynewton.m, так как метод newtons — это просто приложение для итерации с фиксированной точкой, оно построено поверх этого.

person Nicky Mattsson    schedule 02.12.2014
comment
Это не работает. Показывает мне эту ошибку: Error using sym/subsref (line 9) Error using maplemex Error, (in MTM:-subsref) indices must be positive integers, got 229396153.846153855 Error in mynewton>@(x)x-f(x)/fprime(x) (line 8) g = @(x) x - f(x)/fprime(x); Error in fpisimple (line 9) x=f(x); Error in mynewton (line 11) x = fpisimple(g,x0,tol); Error in main (line 51) mynewton(equation,H0,100) - person India Slaver; 02.12.2014