小白快飞 2021-09-08 17:23 采纳率: 66.7%
浏览 108
已结题

【matlab中的subs函数】subs函数运行相关问题

以下是提示的问题:

img


以下是包含与问题相关的代码行:

img

img

img

  • 写回答

1条回答 默认 最新

  • slandarer MATLAB领域优质创作者 2021-09-09 10:58
    关注

    看了一下那本书提供的源代码,其中:
    minRosen.m在第50,80,105行左右第三个参数应用其转置,代码应改为:

    function [x,minf] = minRosen(f,A,b,x0,var,eps)
    format long;
    if nargin == 5
        eps = 1.0e-6;
    end
    
    syms l;
    x0 = transpose(x0);
    n = length(var);
    sz = size(A);
    m = sz(1);
    
    gf = jacobian(f,var);
    bConti = 1;
    
    while bConti
        k = 0;
        s = 0;
        A1 = A;
        A2 = A;
        b1 = b;
        b2 = b;
        for i=1:m
            dfun = A(i,:)*x0 - b(i);
            if abs(dfun)<0.000000001
                k = k + 1;
                A1(k,:) = A(i,:);
                b1(k,1) = b(i);
            else
                s = s+1;
                A2(s,:) = A(i,:);
                b2(s,1) = b(i);
            end
        end
        if k > 0
            A1 = A1(1:k,:);
            b1 = b1(1:k,:);
        end
        if s > 0
            A2 = A2(1:s,:);
            b2 = b2(1:s,:);
        end
        
        while 1
            P = eye(n,n);
            if k > 0
                tM = transpose(A1);
                P = P - tM*inv(A1*tM)*A1;
            end
            gv = Funval(gf, var, x0');
            gv = transpose(gv);
            d = -P*gv;
            if d == 0
                if k == 0
                    x = x0;
                    bConti = 0;
                    break;
                else
                    w = inv(A1*tM)*A1*gv;
                    if w>=0
                        x = x0;
                        bConti = 0;
                        break;
                    else
                        [u,index] = min(w);
                        sA1 = size(A1);
                        if sA1(1) == 1
                            k = 0;
                        else
                            k = sA1(2);
                            A1 =[ A1(1:(index-1),:); A1((index+1):sA1(2),:)];
                        end
                    end
                end
            else
                break;
            end
        end
    
        yl = x0 + l*d;
        tmpf = Funval(f,var,yl');
        bb = b2 - A2*x0;
        dd = A2*d;
        if dd >= 0
            [tmpI,lm] = minJT(tmpf,0,0.1);
        else
            lm = inf;
            for i=1:length(dd)
                if dd(i) < 0
                    if bb(i)/dd(i) < lm
                        lm = bb(i)/dd(i);
                    end
                end
            end
        end
        [xm,minf] = minHJ(tmpf,0,lm,1.0e-14);
        tol = norm(xm*d);
        if tol < eps 
            x = x0;
            break;
        end
        x0 = x0 + xm*d;  
    end
    
    minf = Funval(f,var,x');
    x=double(x);
    minf=double(minf);
    end
    

    为了保证运行速度将minHJ函数x变量取double值:

    function [x,minf] = minHJ(f,a,b,eps)
    format long;
    if nargin == 3
        eps = 1.0e-6;
    end
    
    l = a + 0.382*(b-a);
    u = a + 0.618*(b-a);
    k=1;
    tol = b-a;
    
    while tol>eps && k<100000
        fl = subs(f ,symvar(f), l);
        fu = subs(f ,symvar(f), u);
        if fl > fu
            a = l;
            l = u;
            u = a + 0.618*(b - a);
        else
            b = u;
            u = l;
            l = a + 0.382*(b-a);
        end
        k = k+1;
        tol = abs(b - a);
    end
    if k == 100000
        disp('找不到最小值!');
        x = NaN;
        minf = NaN;
        return;
    end
    x = (a+b)/2;
    x=double(x);
    minf = subs(f, symvar(f),x);
    format short;
    
    

    并将大部分代码中已被移除的函数findsym更改为symvar即可

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论

报告相同问题?

问题事件

  • 系统已结题 9月17日
  • 已采纳回答 9月9日
  • 创建了问题 9月8日

悬赏问题

  • ¥15 求lingo代码和思路
  • ¥15 公交车和无人机协同运输
  • ¥15 stm32代码移植没反应
  • ¥15 matlab基于pde算法图像修复,为什么只能对示例图像有效
  • ¥100 连续两帧图像高速减法
  • ¥15 如何绘制动力学系统的相图
  • ¥15 对接wps接口实现获取元数据
  • ¥20 给自己本科IT专业毕业的妹m找个实习工作
  • ¥15 用友U8:向一个无法连接的网络尝试了一个套接字操作,如何解决?
  • ¥30 我的代码按理说完成了模型的搭建、训练、验证测试等工作(标签-网络|关键词-变化检测)