24小时热门版块排行榜    

查看: 587  |  回复: 0

十甫寸

木虫 (正式写手)

[求助] 拟牛顿法求解非线性方程组——老是出错,跪求大神指点~(或提供一种新的算法)

function [r,m]=mulVNewton(F,x0,A,eps)
%方程组:F
%方程组的初始解:x0
% 初始A矩阵:A
%解的精度:eps
%求得的一组解:r
%迭代步数:m



clear clc
syms x1 x2 x3 f1 f2 f3
D=input('请输入频带的实测数据[q,r,v]:');
q=D(:,1);%第一列
r=D(:,2);
v=D(:,3);

%将实测数据代入方程组中
nn=length(D);
F1=0;F2=0;F3=0;
for i=1:nn
    f1=(v(i)-(x1)*q(i)^(x2)*r(i)^(x3))*q(i)^(x2)*r(i)^(x3);
    f2=(v(i)-(x1)*q(i)^(x2)*r(i)^(x3))*q(i)^(x2)*r(i)^(x3)*log(q(i));
    f3=(v(i)-(x1)*q(i)^(x2)*r(i)^(x3))*q(i)^(x2)*r(i)^(x3)*log(r(i));
    F1=F1+f1;F2=F2+f2;F3=F3+f3;
end
disp('F1=');
disp(F1);
disp('F2=');
disp(F2);
disp('F3=');
disp(F3);

F(1)=F1;F(2)=F2;F(3)=F3;
x=[x1,x2,x3];
x0=[50,0.6,-1.4];
A=[1,0,0; 0,1,0;0,0,1];

if nargin==1
    A=eye(length(x0)); %A取为单位阵
    eps=1.0e-4;
else
   if nargin==2
         eps=1.0e-4;
   end
end

x0 = transpose(x0);
Fx = subs(F, findsym(F),x0);
r=x0-Fx/A;
m=1;
tol=1;
while tol>eps
    x0=r;
    Fx = subs(F, findsym(F),x0);
    r=x0-Fx/A;
    y=r-x0;
    Fr = subs(F, findsym(F),r);
    z= Fr-Fx;
    A1=A+(z-A*y)*transpose(y)/norm(y);                   %调整A
    A=A1;
    m=m+1;
    if(m>100000)                                   %迭代步数控制
        disp('迭代步数太多,可能不收敛!');
        return;
    end
    tol=norm(r-x0);
end

输入数据为[q,r,v]:

[5.6,15.62,1.431
5.6,22.36,1.067
5.6,29.73,0.582
5.6,36.4,0.516
5.6,43.17,0.348
9.4,20.59,1.433
9.4,27.86,0.964
9.4,35.44,0.632
9.4,49.03,0.376
12.8,27.31,1.204
5.6,18.03,1.586
5.6,32.57,0.633
5.6,39.29,0.601
5.6,46.1,0.418
11.6,22.36,1.332
11.6,29.73,0.911
11.6,37.36,0.549
11.6,44.15,0.479
11.6,50.99,0.386]
回复此楼
奋斗改变命运!
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

相关版块跳转 我要订阅楼主 十甫寸 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 没有任何消息-是不是就凉了 +6 图啦图啦 2026-08-24 6/300 2026-08-24 17:44 by 小小红帆
[基金申请] 如果此刻你正在为国基感到焦虑,不妨来听听这首《基金之外》 +5 scalable 2026-08-24 5/250 2026-08-24 17:03 by spiderzhuzhi
[基金申请] 2026国自然函评费到账 +17 羊腰板 2026-08-21 19/950 2026-08-24 16:52 by iaeyuan
[基金申请] 范进中举一文的中心思想 +6 炎黄贵胄 2026-08-22 7/350 2026-08-24 11:58 by 6543yes
[基金申请] 朋友圈看到的 +7 wangzilk 2026-08-18 9/450 2026-08-24 10:50 by cmrandy
[基金申请] 放榜前的不淡定 20+4 snowwithsea 2026-08-19 17/850 2026-08-24 10:20 by echo8914667
[基金申请] 明天应该可查了!? +4 chengyan1220 2026-08-23 4/200 2026-08-24 08:47 by gloomy6159
[基金申请] 2026年的国家社科基金项目通讯评审的新规则与新动向、新挑战 +4 process2012 2026-08-23 5/250 2026-08-23 19:58 by jurkat.1640
[教师之家] 跳槽后在研项目怎么办? +5 简单化xn 2026-08-22 10/500 2026-08-23 12:38 by 简单化xn
[基金申请] 今天放榜吗? +15 布布和一二 2026-08-19 16/800 2026-08-23 09:55 by 张春生
[基金申请] 93BebMhtakh前后11位开头都是大写 +7 且听虎啸 2026-08-17 8/400 2026-08-22 21:55 by 医学老男孩
[基金申请] filecode,4个jtjc了 +13 ziyangfang 2026-08-19 16/800 2026-08-22 17:08 by WH3796
[基金申请] 今天基金会出结果吗?20260819 +16 kkkl_v 2026-08-19 17/850 2026-08-22 16:12 by 阿布Abu
[基金申请] 时间戳今天,20号变了 +5 archvillain 2026-08-20 5/250 2026-08-22 06:12 by hui_daxiao
[基金申请] 看来今天不会放榜了? +8 chengyan1220 2026-08-21 11/550 2026-08-21 17:52 by dcqxinyang
[基金申请] 时间戳又变了 +13 wuchongjun 2026-08-20 19/950 2026-08-21 17:21 by 紫杉醇
[论文投稿] 投稿咨询 +5 wwm09 2026-08-17 7/350 2026-08-21 10:11 by 期刊论文帮手
[基金申请] 重要消息,中午系统在维护 +11 yuleib84 2026-08-18 12/600 2026-08-20 11:09 by xskun
[基金申请] 明天放榜? +5 Shxjjxjkx 2026-08-18 5/250 2026-08-18 18:14 by -大大大大大-
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
信息提示
请填处理意见