matlab根号下是多项式,多项式求根

本帖最后由 pengcsu 于 2016-6-12 12:45 编辑

小弟最近做实验过程中，得到x，y两列数据，数据处理需要得到整数的y对应的x值。但是实验实际测得的都是非整数的y。小弟的思路是对x，y做多项式拟合，得到y(x)表达式，然后求解不同的yn(n=1,2,3...)下，多项式y(x)-yn=0的根，获得整数的yn对应的xn的值。

但是实际编程计算的时候，发现y(x)多项式次数竟然高达20多次，这样一来，迭代初值的选取对于解就很关键了。

于是小弟就编了一个循环语句，从实验数据x，y中找出离yn最接近的y值，用此y对应的x值作为迭代初值来解上述多项式y(x)-yn=0,但此时又提示方程无解。请问高人。上述问题该如何求解？

function chengxu()

clear;clc;

[x]=[24.968

29.968

34.968

39.968

44.968

49.968

54.968

59.968

64.968

69.968

74.968

79.968

84.968

89.968

94.968

99.968

104.968

109.968

114.968

119.968

124.968

129.968

134.968

139.968

144.968

149.968

154.968

159.968

164.968

169.968

174.968

179.968

184.968

189.968

194.968

199.968

204.968

209.968

214.968

219.968

224.968

229.968

234.968

239.968

244.968

249.968

254.968

259.968

264.968

269.968

274.968

279.968

284.968

289.968

294.968

299.968

304.968

309.968

314.968

319.968

324.968

329.968

334.968

339.968

344.968

349.968

354.968

359.968

364.968

369.968

374.968

379.968

384.968

389.968

394.968

399.968

404.968

409.968

414.968

419.968

424.968

429.968

434.968

439.968

444.968

449.968

454.968

459.968

464.968

469.968

474.968

479.968

484.968

489.968

494.968

499.968

504.968

509.968

514.968

519.968

524.968

529.968

534.968

539.968

544.968

549.968

554.968

559.968

564.968

569.968

574.968

579.968

584.968

589.968

594.968

599.968

604.968

609.968

614.968

619.968

624.968

629.968

634.968

639.968

644.968

649.968

654.968

659.968

664.968

669.968

674.968

679.968

684.968

689.968

694.968

699.968

704.968

709.968

714.968

719.968

724.968

729.968

734.968

739.968

744.968

749.968

754.968

759.968

764.968

769.968

774.968

779.968

784.968

789.968

794.968

799.968

804.968

809.968

814.968

819.968

824.968

829.968

834.968

839.968

844.968

849.968

854.968

859.968

864.968

869.968

874.968

879.968

884.968

889.968

894.968

899.968

904.968

909.968

914.968

919.968

924.968

929.968

934.968

939.968

944.968

949.968

954.968

959.968

964.968

969.968

974.968

979.968

984.968

989.968

994.968

]

[y]=[0

0.05472

0.2531

0.36114

0.39585

0.42126

0.4101

0.40093

0.36866

0.34101

0.32258

0.30415

0.29032

0.26055

0.23502

0.22581

0.21198

0.21198

0.22118

0.2212

0.19355

0.18433

0.18433

0.17051

0.15669

0.14286

0.13825

0.12903

0.11521

0.1106

0.10138

0.07373

0.05069

0.04147

0.00921

0.00922

-0.01382

-0.01383

-0.03226

-0.04608

-0.06452

-0.08756

-0.07834

-0.08756

-0.1106

-0.11982

-0.13825

-0.13825

-0.14747

-0.17972

-0.18433

-0.19355

-0.21198

-0.23041

-0.23963

-0.23502

-0.24424

-0.24885

-0.25806

-0.25806

-0.25346

-0.26267

-0.26267

-0.26267

-0.26267

-0.26267

-0.26728

-0.26728

-0.26728

-0.26728

-0.2765

-0.28571

-0.28111

-0.29032

-0.28571

-0.29954

-0.29032

-0.28111

-0.2765

-0.25806

-0.26728

-0.25346

-0.25806

-0.25806

-0.26267

-0.26267

-0.24885

-0.23963

-0.25806

-0.23962

-0.2212

-0.21198

-0.20276

-0.19813

-0.18432

-0.17512

-0.1659

-0.14286

-0.12903

-0.11521

-0.09677

-0.08756

-0.0553

-0.03687

0.01382

0.03226

0.05991

0.09217

0.12903

0.17051

0.19354

0.24419

0.28571

0.32719

0.36871

0.42396

0.47926

0.54377

0.59908

0.66359

0.7235

0.81567

0.8894

0.96322

1.04772

1.12409

1.2212

1.30414

1.39631

1.4779

1.56222

1.67281

1.77875

1.88018

2.00889

2.12442

2.23502

2.341

2.45161

2.5847

2.7098

2.82949

2.9585

3.08756

3.23502

3.38249

3.52535

3.67725

3.82419

3.97678

4.13378

4.29493

4.47467

4.65437

4.82995

5.01374

5.20265

5.40514

5.63134

5.84319

6.08746

6.33656

6.59907

6.89275

7.2212

7.57655

7.95015

8.36406

8.80655

9.29196

9.78712

10.33527

10.91511

11.54455

12.2

12.87039

13.55301

14.27623

14.99236

15.75269

16.5345

17.34496

18.19364

19.04598

19.93502

20.84707

21.79669

22.75194

23.76727

24.74717

25.64362

26.56645

27.54137

28.50983

29.48142

]

m=[0:1:30];

%找出与m值最接近的x，作为迭代初值

for ii=1:length(m)

a=m(ii)

for i=1:195

b(i)=abs(y(i)-a)

[B,n]=min(b)

end

q(ii)=x(n);

end

r=zeros(101,1);

z=zeros(101,101);

%找出最佳拟合次数n

for n=1:100

p=polyfit(x,y,n)

z=polyval(p,x)

r(n)=XGXS(z,y)

[R,n]=max(r);

end

p=polyfit(x,y,n)

hold on

disp('biaozhun')

disp(R)

disp(n)

disp(p)

z=polyval(p,x)

plot(x,y,'ro',x,z);

%% 求解T

T=zeros(length(m),1);

s=zeros(length(m),1);

f=@(T) polyval(p,T);

for ii=1:length(m)

T(ii)=fsolve(@(T) f(T)-m(ii),q(ii));

options=optimset('display','iter','tolfun',1e-2);

end

%% 输出

format long

disp('r为：')

disp(q)

format long

disp(T)

end

%计算相关系数

function r=XGXS(X,Y)

X1=mean(X);

Y1=mean(Y);

fenzi=sum((X-X1).*(Y-Y1));

fenmu=sqrt(sum((X-X1).^2)*sum((Y-Y1).^2));

r=fenzi/fenmu;  %相关系数

end