范文一:matlab数学实验答案
Using Toolbox Path Cache. Type "help toolbox_path_cache" for more info.
To get started, select "MATLAB Help" from the Help menu.
>>
>> clear
>> a=[1,2]
a =
1 2
>> clear
>> a=magic(3)
a =
8 1 6
3 5 7
4 9 2
>> clear
>> z1=2*sin(85*pi/180)/(1+exp(2))
z1 =
0.2375
>> z1=2*sin(85*pi/180)/(1+exp(e)) ??? Undefined function or variable 'e'.
>> x=[2,1+2i;-0.45,5]
z2=(1/2)*log(x+sqrt(1+x*x))
x =
2.0000 1.0000 + 2.0000i
-0.4500 5.0000
z2 =
0.7114 - 0.0253i 0.8968 + 0.3658i
0.2139 + 0.9343i 1.1541 - 0.0044i
>> a=-3.0:0.1:3.0;
z3=((exp(0.3*a)-exp(-0.3*a))/2).*sin(a+0.3)+log((0.3+a)/2)
z3 =
Columns 1 through 5
0.7388 + 3.1416i 0.7696 + 3.1416i 0.7871 + 3.1416i 0.7913 + 3.1416i 0.7822 + 3.1416i
Columns 6 through 10
0.7602 + 3.1416i 0.7254 + 3.1416i 0.6784 + 3.1416i 0.6196 + 3.1416i 0.5496 + 3.1416i
Columns 11 through 15
0.4688 + 3.1416i 0.3780 + 3.1416i 0.2775 + 3.1416i 0.1680 + 3.1416i 0.0497 + 3.1416i
Columns 16 through 20
-0.0771 + 3.1416i -0.2124 + 3.1416i -0.3566 + 3.1416i -0.5104 + 3.1416i -0.6752 + 3.1416i
Columns 21 through 25
-0.8536 + 3.1416i -1.0497 + 3.1416i -1.2701 + 3.1416i -1.5271 + 3.1416i -1.8436 + 3.1416i
Columns 26 through 30
-2.2727 + 3.1416i -2.9837 + 3.1416i -37.0245 -3.0017 -2.3085
Columns 31 through 35
-1.8971 -1.5978 -1.3575 -1.1531 -0.9723
Columns 36 through 40
-0.8083 -0.6567 -0.5151 -0.3819 -0.2561
Columns 41 through 45
-0.1374 -0.0255 0.0792 0.1766
0.2663
Columns 46 through 50
0.3478 0.4206 0.4841 0.5379
0.5815
Columns 51 through 55
0.6145 0.6366 0.6474 0.6470
0.6351
Columns 56 through 60
0.6119 0.5777 0.5327 0.4774
0.4126
Column 61
0.3388
>> t=0:0.5:2.5;
z4=(t.*t).*(t>=0&t<1)+(t.*t-1).*(t>=1&t<2)+(t.*t-2.*t+1).*(t>=2&t<3)>3)>
z4 =
0 0.2500 0 1.2500 1.0000 2.2500
>> A=[12,34,-4;34,7,87;3,65,7];
B=[1,3,-1;2,0,3;3,-2,7];
A+6*B
A-B+ones(size(A))
A*B
A.*B
A^3
A.^3
A/B
A\B
[A,B]
[A([1,3],:);B^2]
ans =
18 52 -10
46 7 105
21 53 49
ans =
12 32 -2
33 8 85
1 68 1
ans =
68 44 62
309 -72 596
154 -5 241
ans =
12 102 4
68 0 261
9 -130 49
ans =
37226 233824 48604
247370 149188 600766
78688 454142 118820
ans =
1728 39304 -64
39304 343 658503
27 274625 343
ans =
16.4000 -13.6000 7.6000
35.8000 -76.2000 50.2000
67.0000 -134.0000 68.0000
ans =
-0.0313 0.3029 -0.3324
0.0442 -0.0323 0.1063
0.0317 -0.1158 0.1558
ans =
12 34 -4 1 3 -1
34 7 87 2 0 3
3 65 7 3 -2 7
ans =
12 34 -4
3 65 7
4 5 1
11 0 19
20 -5 40
>> A=[1,2,3,4,5;6,7,8,9,10;11,12,13,14,15;16,17,18,19,20;21,22,23,24,25];
B=[3,0,16;17,-6,9;0,23,-4;9,7,0;4,13,11];
C=A*B
D=C(3:end,2:end)
whos
C =
93 150 77
258 335 237
423 520 397
588 705 557
753 890 717
D =
520 397
705 557
890 717
Name Size Bytes Class
A 5x5 200 double array
B 5x3 120 double array
C 5x3 120 double array
D 3x2 48 double array
a 1x61 488 double array
ans 5x3 120 double array
t 1x6 48 double array
x 2x2 64 double array (complex)
z1 1x1 8 double array
z2 2x2 64 double array (complex)
z3 1x61 976 double array (complex)
z4 1x6 48 double array
Grand total is 219 elements using 2304 bytes
>> 最后一题x=100:999;
length(find(rem(x,21)==0))
??? 最后一题x=100:999;
|
Error: Missing variable or function.
>> clear
>> X=100:999;
>> length(find(rem(x,21)==0))
??? Undefined function or variable 'x'.
>> X=100:1:999
X =
Columns 1 through 11
100 101 102 103 104 105 106 107 108 109 110
Columns 12 through 22
111 112 113 114 115 116 117 118 119 120 121
Columns 23 through 33
122 123 124 125 126 127 128 129 130 131 132
Columns 34 through 44
133 134 135 136 137 138 139 140 141 142 143
Columns 45 through 55
144 145 146 147 148 149 150 151 152 153 154
Columns 56 through 66
155 156 157 158 159 160 161 162 163 164 165
Columns 67 through 77
166 167 168 169 170 171 172 173 174 175 176
Columns 78 through 88
177 178 179 180 181 182 183 184 185 186 187
Columns 89 through 99
188 189 190 191 192 193 194 195 196 197 198
Columns 100 through 110
199 200 201 202 203 204 205 206 207 208 209
Columns 111 through 121
210 211 212 213 214 215 216 217 218 219 220
Columns 122 through 132
221 222 223 224 225 226 227 228 229 230 231
Columns 133 through 143
232 233 234 235 236 237 238 239 240 241 242
Columns 144 through 154
243 244 245 246 247 248 249 250 251 252 253
Columns 155 through 165
254 255 256 257 258 259 260 261 262 263 264
Columns 166 through 176
265 266 267 268 269 270 271 272 273 274 275
Columns 177 through 187
276 277 278 279 280 281 282 283 284 285 286
Columns 188 through 198
287 288 289 290 291 292 293 294 295 296 297
Columns 199 through 209
298 299 300 301 302 303 304 305 306 307 308
Columns 210 through 220
309 310 311 312 313 314 315 316 317 318 319
Columns 221 through 231
320 321 322 323 324 325 326 327 328 329 330
Columns 232 through 242
331 332 333 334 335 336 337 338 339 340 341
Columns 243 through 253
342 343 344 345 346 347 348 349 350 351 352
Columns 254 through 264
353 354 355 356 357 358 359 360 361 362 363
Columns 265 through 275
364 365 366 367 368 369 370 371 372 373 374
Columns 276 through 286
375 376 377 378 379 380 381 382 383 384 385
Columns 287 through 297
386 387 388 389 390 391 392 393 394 395 396
Columns 298 through 308
397 398 399 400 401 402 403 404 405 406 407
Columns 309 through 319
408 409 410 411 412 413 414 415 416 417 418
Columns 320 through 330
419 420 421 422 423 424 425 426 427 428 429
Columns 331 through 341
430 431 432 433 434 435 436 437 438 439 440
Columns 342 through 352
441 442 443 444 445 446 447 448 449 450 451
Columns 353 through 363
452 453 454 455 456 457 458 459 460 461 462
Columns 364 through 374
463 464 465 466 467 468 469 470 471 472 473
Columns 375 through 385
474 475 476 477 478 479 480 481 482 483 484
Columns 386 through 396
485 486 487 488 489 490 491 492 493 494 495
Columns 397 through 407
496 497 498 499 500 501 502 503 504 505 506
Columns 408 through 418
507 508 509 510 511 512 513 514 515 516 517
Columns 419 through 429
518 519 520 521 522 523 524 525 526 527 528
Columns 430 through 440
529 530 531 532 533 534 535 536 537 538 539
Columns 441 through 451
540 541 542 543 544 545 546 547 548 549 550
Columns 452 through 462
551 552 553 554 555 556 557 558 559 560 561
Columns 463 through 473
562 563 564 565 566 567 568 569 570 571 572
Columns 474 through 484
573 574 575 576 577 578 579 580 581 582 583
Columns 485 through 495
584 585 586 587 588 589 590 591 592 593 594
Columns 496 through 506
595 596 597 598 599 600 601 602 603 604 605
Columns 507 through 517
606 607 608 609 610 611 612 613 614 615 616
Columns 518 through 528
617 618 619 620 621 622 623 624 625 626 627
Columns 529 through 539
628 629 630 631 632 633 634 635 636 637 638
Columns 540 through 550
639 640 641 642 643 644 645 646 647 648 649
Columns 551 through 561
650 651 652 653 654 655 656 657 658 659 660
Columns 562 through 572
661 662 663 664 665 666 667 668 669 670 671
Columns 573 through 583
672 673 674 675 676 677 678 679 680 681 682
Columns 584 through 594
683 684 685 686 687 688 689 690 691 692 693
Columns 595 through 605
694 695 696 697 698 699 700 701 702 703 704
Columns 606 through 616
705 706 707 708 709 710 711 712 713 714 715
Columns 617 through 627
716 717 718 719 720 721 722 723 724 725 726
Columns 628 through 638
727 728 729 730 731 732 733 734 735 736 737
Columns 639 through 649
738 739 740 741 742 743 744 745 746 747 748
Columns 650 through 660
749 750 751 752 753 754 755 756 757 758 759
Columns 661 through 671
760 761 762 763 764 765 766 767 768 769 770
Columns 672 through 682
771 772 773 774 775 776 777 778 779 780 781
Columns 683 through 693
782 783 784 785 786 787 788 789 790 791 792
Columns 694 through 704
793 794 795 796 797 798 799 800 801 802 803
Columns 705 through 715
804 805 806 807 808 809 810 811 812 813 814
Columns 716 through 726
815 816 817 818 819 820 821 822 823 824 825
Columns 727 through 737
826 827 828 829 830 831 832 833 834 835 836
Columns 738 through 748
837 838 839 840 841 842 843 844 845 846 847
Columns 749 through 759
848 849 850 851 852 853 854 855 856 857 858
Columns 760 through 770
859 860 861 862 863 864 865 866 867 868 869
Columns 771 through 781
870 871 872 873 874 875 876 877 878 879 880
Columns 782 through 792
881 882 883 884 885 886 887 888 889 890 891
Columns 793 through 803
892 893 894 895 896 897 898 899 900 901 902
Columns 804 through 814
903 904 905 906 907 908 909 910 911 912 913
Columns 815 through 825
914 915 916 917 918 919 920 921 922 923 924
Columns 826 through 836
925 926 927 928 929 930 931 932 933 934 935
Columns 837 through 847
936 937 938 939 940 941 942 943 944 945 946
Columns 848 through 858
947 948 949 950 951 952 953 954 955 956 957
Columns 859 through 869
958 959 960 961 962 963 964 965 966 967 968
Columns 870 through 880
969 970 971 972 973 974 975 976 977 978 979
Columns 881 through 891
980 981 982 983 984 985 986 987 988 989 990
Columns 892 through 900
991 992 993 994 995 996 997 998 999
>> x=100:1:999;
length(find(rem(x,21)==0))
ans =
43
>>
范文二:matlab数学实验答案
%Page20,ex1
(5) 等于 [exp(1),exp(2);exp(3),exp(4)]
(7) 3=1*3, 8=2*4
(8) a为各列最小值, b 为最小值所在的行号
(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture
(11) 答案表明:编址第 2元素满足不等式 (30>=20)和编址第 4元素满足不等式 (40>=10) (12) 答案表明:编址第 2行第 1列元素满足不等式 (30>=20)和编址第 2行第 2列元素满足不 等式 (40>=10)
%Page20, ex2
(1)a, b, c的值尽管都是 1, 但数据类型分别为数值,字符, 逻辑, 注意 a 与 c 相等, 但 他们不等于 b
(2)double(fun)输出的分别是字符 a,b,s,(,x,)的 ASCII 码
%Page20,ex3
>> r=2;p=0.5;n=12;
>> T=log(r)/n/log(1+0.01*p)
T =
11.5813
%Page20,ex4
>> x=-2:0.05:2;f=x.^4-2.^x;
>> [fmin,min_index]=min(f)
fmin =
-1.3907 %最小值
min_index =
54 %最小值点编址
>> x(min_index)
ans =
0.6500 %最小值点
>> [f1,x1_index]=min(abs(f)) %求近似根 --绝对值最小的点
f1 =
0.0328
x1_index =
24
>> x(x1_index)
ans =
-0.8500
>> x(x1_index)=[];f=x.^4-2.^x; %删去绝对值最小的点以求函数绝对值次小的点
>> [f2,x2_index]=min(abs(f)) %求另一近似根 --函数绝对值次小的点
f2 =
0.0630
x2_index =
65
>> x(x2_index)
ans =
1.2500
%Page20,ex5
>> z=magic(10)
z =
92 99 1 8 15 67 74 51 58 40
98 80 7 14 16 73 55 57 64 41
4 81 88 20 22 54 56 63 70 47
85 87 19 21 3 60 62 69 71 28
86 93 25 2 9 61 68 75 52 34
17 24 76 83 90 42 49 26 33 65
23 5 82 89 91 48 30 32 39 66
79 6 13 95 97 29 31 38 45 72
10 12 94 96 78 35 37 44 46 53
11 18 100 77 84 36 43 50 27 59
>> sum(z)
ans =
505 505 505 505 505 505 505 505 505 505 >> sum(diag(z))
ans =
505
>> z(:,2)/sqrt(3)
ans =
57.1577
46.1880
46.7654
50.2295
53.6936
13.8564
2.8868
3.4641
6.9282
10.3923
>> z(8,:)=z(8,:)+z(3,:)
z =
92 99 1 8 15 67 74 51 58 40
98 80 7 14 16 73 55 57 64 41
4 81 88 20 22 54 56 63 70 47
85 87 19 21 3 60 62 69 71 28
86 93 25 2 9 61 68 75 52 34
17 24 76 83 90 42 49 26 33 65
23 5 82 89 91 48 30 32 39 66
83 87 101 115 119 83 87 101 115 119 10 12 94 96 78 35 37 44 46 53
11 18 100 77 84 36 43 50 27 59
%Page 40 ex1
先在编辑器窗口写下列 M 函数,保存为 eg2_1.m function [xbar,s]=ex2_1(x)
n=length(x);
xbar=sum(x)/n;
s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));
例如
>>x=[81 70 65 51 76 66 90 87 61 77];
>>[xbar,s]=ex2_1(x)
xbar =
72.4000
s =
12.1124
%Page 40 ex2
s=log(1);n=0;
while s<>
n=n+1;
s=s+log(1+n);
end
m=n
计算结果 m=37
%Page 40 ex3
clear;
F(1)=1;F(2)=1;k=2;x=0;
e=1e-8; a=(1+sqrt(5))/2;
while abs(x-a)>e
k=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);
end
a,x,k
计算至 k=21可满足精度
%Page 40 ex4
clear;tic;s=0;
for i=1:1000000
s=s+sqrt(3)/2^i;
end
s,toc
tic;s=0;i=1;
while i<>
s=s+sqrt(3)/2^i;i=i+1;
end
s,toc
tic;s=0;
i=1:1000000;
s=sqrt(3)*sum(1./2.^i);
s,toc
%Page 40 ex5
t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ... 31 32 31 29 27 25 24 22 20 18 17 16];
plot(t,c)
%Page 40 ex6
%(1)
x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y) y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2]) %(2)参数方法
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t); plot(x,y)
%(3)
x=-3:0.1:3;y=x;
[x,y]=meshgrid(x,y);
z=x.^2+y.^2;
surf(x,y,z)
%(4)
x=-3:0.1:3;y=-3:0.1:13;
[x,y]=meshgrid(x,y);
z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6; surf(x,y,z)
%(5)
t=0:0.01:2*pi;
x=sin(t);y=cos(t);z=cos(2*t);
plot3(x,y,z)
%(6)
theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20); [theta,fai]=meshgrid(theta,fai);
x=2*sin(fai).*cos(theta);
y=2*sin(fai).*sin(theta);z=2*cos(fai);
surf(x,y,z)
%(7)
x=linspace(0,pi,100);
y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);
plot(x,y1,x,y2,x,y3)
%page41, ex7
x=-1.5:0.05:1.5;
y=1.1*(x>1.1)+x.*(x<=1.1).*(x>=-1.1)-1.1*(x<>
plot(x,y)
%page41,ex8
分别使用 which trapz, type trapz, dir C:\MATLAB7\toolbox\matlab\datafun\ %page41,ex9
clear;close;
x=-2:0.1:2;y=x;
[x,y]=meshgrid(x,y);
a=0.5457;b=0.7575;
p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);
p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y<>
p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<>
mesh(x,y,p)
%page41, ex10
lookfor lyapunov
help lyap
>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];
>> X=lyap(A,C)
X =
1.0000 -1.0000 -0.0000
-1.0000 2.0000 1.0000
-0.0000 1.0000 7.0000
%Chapter 3
%Exercise 1
>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\b
ans =
0.5000 0.5000 1.0000
ans =
2 2 1
ans =
0.6552 %一元方程组 x[2,4,3]=[1,2,3]的近似解
ans =
0 0 0
0 0 0
0.6667 1.3333 1.0000
%矩阵方程 [1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解 %Exercise 2(1)
>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];
>> rank(A), rank([A,b]) %[A,b]为增广矩阵
ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
2.3830
1.4894
2.0213
%Exercise 2(2)
>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];
>> rank(A), rank([A,b])
ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
-0.4706
-0.2941
%Exercise 2(3)
>> A=[4 1;3 2;1 -5];b=[1;1;1];
>> rank(A), rank([A,b])
ans =
2
ans =
3 %可见方程组无解
>> x=A\b
x =
0.3311
-0.1219 %最小二乘近似解
%Exercise 2(4)
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意 b 的写法 >> rank(a),rank([a,b])
ans =
3
ans =
3 %rank(a)==rank([a,b])<>
>> a\b
ans =
1
1
0 %一个特解
%Exercise 3
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';
>> x=null(a),x0=a\b
x =
-0.6255
0.6255
-0.2085
0.4170
x0 =
1
1
%通解 kx+x0
%Exercise 4
>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95];
>> x1=a*x, x2=a^2*x, x10=a^10*x
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> x0=[0.8 0.2]';
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> [v,e]=eig(a)
v =
0.9806 -0.7071
0.1961 0.7071
e =
1.0000 0
0 0.9400
>> v(:,1)./x
ans =
1.1767
1.1767 %成比例,说明 x 是最大特征值对应的特征向量 %Exercise 5
%用到公式 (3.11)(3.12)
>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]';
>> C=B/diag(x)
C =
0.2400 0.4000 0.0500
0.0900 0.2000 0.0100
0.1200 0.0400 0.0900
>> A=eye(3,3)-C
A =
0.7600 -0.4000 -0.0500
-0.0900 0.8000 -0.0100
-0.1200 -0.0400 0.9100
>> D=[17 17 17]';x=A\D
x =
37.5696
25.7862
24.7690
%Exercise 6(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans =
-94
ans =
0.2553 -0.0213 0.0426
0.1596 -0.1383 -0.2234
0.1809 -0.2234 -0.0532
v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
%Exercise 6(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans =
1
ans =
2.0000 -2.0000 1.0000
1.0000 -1.0000 1.0000
2.0000 -3.0000 2.0000
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774
-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0
0 0 1.0000 - 0.0000i
%Exercise 6(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> det(A),inv(A), [v,d]=eig(A)
ans =
1
ans =
68.0000 -41.0000 -17.0000 10.0000
-41.0000 25.0000 10.0000 -6.0000
-17.0000 10.0000 5.0000 -3.0000
10.0000 -6.0000 -3.0000 2.0000
v =
0.8304 0.0933 0.3963 0.3803
-0.5016 -0.3017 0.6149 0.5286
-0.2086 0.7603 -0.2716 0.5520
0.1237 -0.5676 -0.6254 0.5209
d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
%Exercise 6(4)(以 n=5为例 )
%关键是矩阵的定义
%方法一(三个 for )
n=5;
for i=1:n, a(i,i)=5;end
for i=1:(n-1),a(i,i+1)=6;end
for i=1:(n-1),a(i+1,i)=1;end
a
%方法二(一个 for )
n=5;a=zeros(n,n);
a(1,1:2)=[5 6];
for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];end
a(n,[n-1 n])=[1 5];
a
%方法三(不用 for )
n=5;a=diag(5*ones(n,1));
b=diag(6*ones(n-1,1));
c=diag(ones(n-1,1));
a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)] %下列计算
>> det(a)
ans =
665
>> inv(a)
ans =
0.3173 -0.5865 1.0286 -1.6241 1.9489
-0.0977 0.4887 -0.8571 1.3534 -1.6241
0.0286 -0.1429 0.5429 -0.8571 1.0286
-0.0075 0.0376 -0.1429 0.4887 -0.5865
0.0015 -0.0075 0.0286 -0.0977 0.3173
>> [v,d]=eig(a)
v =
-0.7843 -0.7843 -0.9237 0.9860 -0.9237
0.5546 -0.5546 -0.3771 -0.0000 0.3771
-0.2614 -0.2614 0.0000 -0.1643 0.0000
0.0924 -0.0924 0.0628 -0.0000 -0.0628
-0.0218 -0.0218 0.0257 0.0274 0.0257
d =
0.7574 0 0 0 0
0 9.2426 0 0 0
0 0 7.4495 0 0
0 0 0 5.0000 0
0 0 0 0 2.5505
%Exercise 7(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a)
v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
>> det(v)
ans =
-0.9255 %v行列式正常 , 特征向量线性相关,可对角化 >> inv(v)*a*v %验算
ans =
-3.0527 0.0000 -0.0000
0.0000 3.6760 -0.0000
-0.0000 -0.0000 8.3766
>> [v2,d2]=jordan(a) %也可用 jordan
v2 =
0.0798 0.0076 0.9127
0.1886 -0.3141 0.1256
-0.1605 -0.2607 0.4213 %特征向量不同
d2 =
8.3766 0 0
0 -3.0527 - 0.0000i 0
0 0 3.6760 + 0.0000i
>> v2\a*v2
ans =
8.3766 0 0.0000
0.0000 -3.0527 0.0000
0.0000 0.0000 3.6760
>> v(:,1)./v2(:,2) %对应相同特征值的特征向量成比例
ans =
2.4491
2.4491
2.4491
%Exercise 7(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774
-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0
0 0 1.0000 - 0.0000i
>> det(v)
ans =
-5.0566e-028 -5.1918e-017i %v的行列式接近 0, 特征向量线性相关,不可对角化 >> [v,d]=jordan(a)
v =
1 0 1
1 0 0
1 -1 0
d =
1 1 0
0 1 1
0 0 1 %jordan标准形不是对角的,所以不可对角化
%Exercise 7(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> [v,d]=eig(A)
v =
0.8304 0.0933 0.3963 0.3803
-0.5016 -0.3017 0.6149 0.5286
-0.2086 0.7603 -0.2716 0.5520
0.1237 -0.5676 -0.6254 0.5209
d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
>> inv(v)*A*v
ans =
0.0102 0.0000 -0.0000 0.0000
0.0000 0.8431 -0.0000 -0.0000
-0.0000 0.0000 3.8581 -0.0000
-0.0000 -0.0000 0 30.2887
%本题用 jordan 不行 , 原因未知
%Exercise 7(4)参考 6(4)和 7(1), 略
%Exercise 8 只有 (3)对称 , 且特征值全部大于零 , 所以是正定矩阵 . %Exercise 9(1)
>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0]
>> rank(a)
ans =
3
>> rank(a(1:3,:))
ans =
2
>> rank(a([1 2 4],:)) %1,2,4行为最大无关组
ans =
3
>> b=a([1 2 4],:)';c=a([3 5],:)';
>> b\c %线性表示的系数
ans =
0.5000 5.0000
-0.5000 1.0000
0 -5.0000
%Exercise 10
>> a=[1 -2 2;-2 -2 4;2 4 -2]
>> [v,d]=eig(a)
v =
0.3333 0.9339 -0.1293
0.6667 -0.3304 -0.6681
-0.6667 0.1365 -0.7327
d =
-7.0000 0 0
0 2.0000 0
0 0 2.0000
>> v'*v
ans =
1.0000 0.0000 0.0000
0.0000 1.0000 0
0.0000 0 1.0000 %v确实是正交矩阵
%Exercise 11
%设经过 6个电阻的电流分别为 i1, ..., i6. 列方程组如下
%20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;
%i1=i3+i4;i5=i2+i3;i6=i4+i5;
%计算如下
>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0; 0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0; 0 0 0 0 0 0 -1 -1 1];
>>b=[20 5 0 0 0 0 0 0 0]'; A\b
ans =
13.3453
6.4401
8.5420
3.3274
-1.1807
1.6011
1.7263
0.4204
2.1467
%Exercise 12
>> A=[1 2 3;4 5 6;7 8 0];
>> left=sum(eig(A)), right=sum(trace(A))
left =
6.0000
right =
6
>> left=prod(eig(A)), right=det(A) %原题有错 , (-1)^n应删去
left =
27.0000
right =
27
>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3)) fA =
1.0e-012 *
0.0853 0.1421 0.0284
0.1421 0.1421 0
-0.0568 -0.1137 0.1705
>> norm(fA) %f(A)范数接近 0
ans =
2.9536e-013
%Exercise 1(1)
roots([1 1 1])
%Exercise 1(2)
roots([3 0 -4 0 2 -1])
%Exercise 1(3)
p=zeros(1,24);
p([1 17 18 22])=[5 -6 8 -5];
roots(p)
%Exercise 1(4)
p1=[2 3];
p2=conv(p1, p1);
p3=conv(p1, p2);
p3(end)=p3(end)-4; %原 p3最后一个分量 -4
roots(p3)
%Exercise 2
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
fzero(fun,2)
%Exercise 3
fun=inline('x^4-2^x');
fplot(fun,[-2 2]);grid on;
fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)
%Exercise 4
fun=inline('x*sin(1/x)','x');
fplot(fun, [-0.1 0.1]);
x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;
x=[x,-x]
%Exercise 5
fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^ 2-16*x(3)^2]','x');
[a,b,c]=fsolve(fun,[0 0 0])
%Exercise 6
fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];
[a,b,c]=fsolve(fun,[0.5 0.5])
%Exercise 7
clear; close; t=0:pi/100:2*pi;
x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);
x2=3+sqrt(2)*cos(t); y2=6*sin(t);
plot(x1,y1,x2,y2); grid on; %作图发现 4个解的大致位置,然后分别求解
y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])
y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])
y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])
y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])
%Exercise 8(1)
clear;
fun=inline('x.^2.*sin(x.^2-x-2)');
fplot(fun,[-2 2]);grid on; %作图观察
x(1)=-2;
x(3)=fminbnd(fun,-1,-0.5);
x(5)=fminbnd(fun,1,2);
fun2=inline('-x.^2.*sin(x.^2-x-2)');
x(2)=fminbnd(fun2,-2,-1);
x(4)=fminbnd(fun2,-0.5,0.5);
x(6)=2
feval(fun,x)
%答案 : 以上 x(1)(3)(5)是局部极小, x(2)(4)(6)是局部极大,从最后一句知道 x(1)全局最小, x(2)最大。
%Exercise 8(2)
clear;
fun=inline('3*x.^5-20*x.^3+10');
fplot(fun,[-3 3]);grid on;%作图观察
x(1)=-3;
x(3)=fminsearch(fun,2.5);
fun2=inline('-(3*x.^5-20*x.^3+10)');
x(2)=fminsearch(fun2,-2.5);
x(4)=3;
feval(fun,x)
%Exercise 8(3)
fun=inline('abs(x^3-x^2-x-2)');
fplot(fun,[0 3]);grid on;%作图观察
fminbnd(fun,1.5,2.5)
fun2=inline('-abs(x^3-x^2-x-2)');
fminbnd(fun2,0.5,1.5)
%Exercise 9
close;
x=-2:0.1:1;y=-7:0.1:1;
[x,y]=meshgrid(x,y);
z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9;
mesh(x,y,z);grid on;%作图观察
fun=inline('x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9');
x=fminsearch(fun,[0 0])%求极小值
fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)');
x=fminsearch(fun2,[0 -5])%求极大值
%Exercise 10
clear;t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16];
p2=polyfit(t,c,2)
p3=polyfit(t,c,3)
fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t');
a=lsqcurvefit(fun,[0 0],t,c)%初值可以试探
f=feval(fun, a,t)
norm(f-c)%拟合效果
plot(t,c,t,f) %作图检验
fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改 f(x)+20 b=lsqcurvefit(fun2,[0 0],t,c)
figure
f2=feval(fun2, b,t)
norm(f2-c)%拟合效果
plot(t,c,t,f2) %作图检验
%Exercise 11
fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)'); x=fzero(fun, 0, 1)
%Exercise 12
r=5.04/12/100;N=20*12;
x=7500*180 %房屋总价格
y=x*0.3 %首付款额
x0=x-y%贷款总额
a=(1+r)^N*r*x0/((1+r)^N-1)%月付还款额
r1=4.05/12/100;x1=10*10000;%公积金贷款
a1=(1+r1)^N*r1*x1/((1+r1)^N-1)
x2=x0-x1%商业贷款
a2=(1+r)^N*r*x2/((1+r)^N-1)
a=a1+a2
%Exercise 13
%列方程 th*R^2+(pi-2*th)*r^2-R*r*sin(th)=pi*r^2/2
%化简得 sin(2*th)-2*th*cos(2*th)=pi/2
%以下 Matlab 计算
clear;fun= inline('sin(2*th)-2*th*cos(2*th)-pi/2','th')
th=fsolve(fun,pi/4)
R=20*cos(th)
%Exercise 14
%先在 Editor 窗口写 M 函数保存
function x=secant(fname,x0,x1,e)
while abs(x0-x1)>e,
x=x1-(x1-x0)*feval(fname,x1)/(feval(fname,x1)-feval(fname,x0));
x0=x1;x1=x;
end
%再在指令窗口
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x'); secant(fun,1,2,1e-8)
%Exercise 15
%作系数为 a, 初值为 xo, 从第 m 步到第 n 步迭代过程的 M 函数: function f=ex4_15fun(a,x0,m,n)
x(1)=x0; y(1)=a*x(1)+1;x(2)=y(1);
if m<2, plot([x(1),x(1),x(2)],[0,y(1),y(1)]);hold="" on;="">2,>
for i=2:n
y(i)=a*x(i)+1; x(i+1)=y(i);
if i>m, plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]); end
end
hold off;
%M脚本文件
subplot(2,2,1);ex4_15fun(0.9,1,1,20);
subplot(2,2,2);ex4_15fun(-0.9,1,1,20);
subplot(2,2,3);ex4_15fun(1.1,1,1,20);
subplot(2,2,4);ex4_15fun(-1.1,1,1,20);
%Exercise 16
%设夹角 t, 问题转化为 min f=5/sin(t)+10/cos(t)
%取初始值 pi/4, 计算如下
fun=@(t)5/sin(t)+10/cos(t);
[t,f]=fminsearch(fun, pi/4)
t =
0.6709
f =
20.8097
%Exercise 17
%提示:x(k+2)=f(x(k))=a^2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k))) %计算平衡点 x
%|f'(x)|<>
%Exercise 18
%先写 M 文件
function f=ex4_18(a,x0,n)
x=zeros(1,n);y=x;
x(1)=x0;
y(1)=a*x(1)+1;
x(2)=y(1);
plot([x(1),x(1),x(2)],[0,y(1),y(1)],'r');
hold on;
for i=2:n
y(i)=a*x(i)+1;
x(i+1)=y(i);
plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]) end
hold off;
%再执行指令
>> ex4_18(0.9,1,20)
>> ex4_18(-0.9,1,20)
>> ex4_18(1.1,1,20)
>> ex4_18(-1.1,1,20)
%Exercise 19
clear; close; x(1)=0; y(1)=0;
for k=1:3000
x(k+1)=1+y(k)-1.4*x(k)^2; y(k+1)=0.3*x(k); end
plot(x(1000:1500),y(1000:1500),'+g');hold on plot(x(1501:2000),y(1501:2000),'.b'); plot(x(2001:2500),y(2001:2500),'*y'); plot(x(2501:3001),y(2501:3001),'.r');
%Exercise 1
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0];
trapz(x,y)
%Exercise 2
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0];
diff(y)./diff(x)
%Exercise 3
xa=-1:0.1:1;ya=0:0.1:2;
[x,y]=meshgrid(xa,ya);
z=x.*exp(-x.^2 -y.^3);
[px,py] = gradient(z,xa,ya);
px
%Exercise 4
t=0:0.01:1.5;
x=log(cos(t));
y=cos(t)-t.*sin(t);
dydx=gradient(y,x)
[x_1,id]=min(abs(x-(-1)));%找最接近 x=-1的点 dydx(id)
%Exercise 5(2)
fun=inline('exp(2*x).*cos(x).^3'); quadl(fun,0,2*pi)
或用 trapz
x=linspace(0,2*pi,100);
y=exp(2*x).*cos(x).^3;
trapz(x,y)
%Exercise 5(3)
fun=@(x)x.*log(x.^4).*asin(1./x.^2); quadl(fun,1,3)
或用 trapz
x=1:0.01:3;
y=feval(fun,x);
trapz(x,y)
%Exercise 5(4)
fun=@(x)sin(x)./x;
quadl(fun,1e-10,1) %注意由于下限为 0,被积函数没有意义,用很小的 1e-10代替 %Exercise 5(5)
%参考 Exercise 5(4)
%Exercise 5(6)
fun=inline('sqrt(1+r.^2.*sin(th))','r','th');
dblquad(fun,0,1,0,2*pi)
%Exercise 5(7)
首先建立 84页函数 dblquad2
clear;
fun=@(x,y)1+x+y.^2;
clo=@(x)-sqrt(2*x-x.^2);
dup=@(x)sqrt(2*x-x.^2);
dblquad2(fun,0,2,clo,dhi,100)
%Exercise 6
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t);
dx=gradient(x,t);dy=gradient(y,t);
f=sqrt(dx.^2+dy.^2);
trapz(t,f)
%Exercise 7
xa=-1:0.1:1;ya=0:0.1:2;
[x,y]=meshgrid(xa,ya);
z=x.*exp(x.^2+y.^2);
[zx,zy]=gradient(z,xa,ya);
f=sqrt(1+zx.^2+zy.^2);
s=0;
for i=2:length(xa)
for j=2:length(ya)
s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4;
end
end
s
%Exercise 8
funl=inline('-(-x).^0.2.*cos(x)');
funr=inline('x.^0.2.*cos(x)');
quadl(funl,-1,0)+quadl(funr,0,1)
%Exercise 9 (以 I32为例 )
fun=@(x)abs(sin(x));
h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y)
h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致,结果失真 q1=quad(fun,0,32*pi)
q2=quadl(fun,0,32*pi)
%Exercise 10(2)
%先在程序编辑器,写下列函数,保存为 ex5_10_2f
function d=ex5_10_2f(fname,a,h0,e)
h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h); d0=d+2*e;
while abs(d-d0)>e
d0=d;h0=h;h=h0/2;
d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);
end
%再在指令窗口执行
fun=inline('x.^2*sin(x.^2-x-2)','x');
d=ex5_10_2f(fun,1.4,0.1,1e-3)
%Exercise 11
%提示:f 上升时, f'>0;f下降时, f'<0; f极值,="" f'="">0;>
%Exercise 12
在程序编辑器,写下列函数,保存为 ex5_12f
function I=ex5_12(fname,a,b,n)
x=linspace(a,b,n+1);
y=feval(fname,x);
I=(b-a)/n/3*(y(1)+y(n+1)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));
%再在指令窗口执行
ex5_12(inline('1/sqrt(2*pi)*exp(-x.^2/2)'),0,1,50)
%Exercise 13
fun=inline('5400*v./(8.276*v.^2+2000)','v');
quadl(fun,15,30)
%Exercise 14
重心不超过凳边沿。 1/2, 2/3, 3/4, ...,n/(n+1)
%Exercise15
利润函数 fun=inline('(p-c0+k*log(M*exp(-a*p)))*M*exp(-a*p)','p'); 求 p 使 fun 最大
%Exercise 16
clear; x=-3/4:0.01:3/4;
y=(3/4+x)*2.*sqrt(1-16/9.*x.^2)*9.8;
P=trapz(x,y) %单位:千牛
%Exercise 17
clear; close;
fplot('17-t^(2/3)-5-2*t^(2/3)',[0,20]); grid;
t=fzero('17-x^(2/3)-5-2*x^(2/3)',7)
t=0:0.1:8; y=17-t.^(2/3)-5-2*t.^(2/3);
trapz(t,y)-20 %单位:百万元
%Exercise 18
%曲面面积计算
%Excercise 1(1)
fun=inline('x+y','x','y');
[t,y]=ode45(fun,[0 1 2 3],1) %注意由于初值为 y(0)=1,[0 1 2 3]中 0不可缺
%Excercise 1(3)
%令 y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=0.01*y(2)^2-2*y(1)+sin(t)
%运行下列指令
clear;close;
fun=@(t,y)[y(2);0.01*y(2)^2-2*y(1)+sin(t)];
[t,y]=ode45(fun,[0 5],[0;1]);
plot(t,y(:,1))
%Excercise 1(5)
%令 y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=-mu*(y(1)^2-1)*y(2)-y(1)
%运行下列指令 , 注意参数 mu 的处理
clear;close;
fun=@(t,y,mu)[y(2);-mu*(y(1)^2-1)*y(2)-y(1)];
[t,y]=ode45(fun,[0 20],[2;0],[],1);
plot(y(:,1),y(:,2));hold on;
[t,y]=ode45(fun,[0 20],[2;0],[],2);
plot(y(:,1),y(:,2),'r');hold off;
%Excercise 2
roots([1 10 54 132 137 50])
%通解 A1*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-t)+A5*t*exp(-t) %Excercise 3
dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');
[x,y]=ode45(dfun,[0,50],[1;-1]);length(x)
%所用节点很多
[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x)
%所用节点很少
%Excercise 4
clear;
dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-( 1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3;
x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1 )+1-1/82.45)^2+x(3)^2))^3]','t','x');
[t,x]=ode45(dfun,[0 24],[1.2; 0; 0; -1.04935371]);
plot(x(:,1),x(:,3));
%Excercise 5
%方程 y'=2x+y^2,y(0)=0
clear;close;
fun=inline('2*x+y^2','x','y');
[x,y]=ode45(fun,[0 1.57],0); %x的上界再增加 , 解会
plot(x,y)
%Excercise 6
clear;close;
fun=@(t,x,a,b)a*x+b;
[t,x]=ode45(fun,[0 10],0.1,[],1,1); subplot(2,4,1);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,1); subplot(2,4,2);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],1,-1); subplot(2,4,3);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,-1); subplot(2,4,4);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,1); subplot(2,4,5);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,1); subplot(2,4,6);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,-1); subplot(2,4,7);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,-1); subplot(2,4,8);plot(t,x) %Excercise 7
%微分方程 T'=k(c-T),T(0)=20 dsolve('DT=k*(c-T)','T(0)=20','t') %得 c+exp(-k*t)*(-c+20)
%利用 T(10)=25.2, T(20)=28.32拟合 (或者解非线性方程 ) fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t') lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32])
%解得户外温度 c=33,比例系数 k=0.05.
%Excercise 8
%微分方程 x'/x=0.5*(1-x),x(0)=0.05
fun=inline('0.5*(1-x)*x','t','x');
[t,x]=ode45(fun,[0 10],0.05);
plot(t,x)
id=min(find(x>0.8));
t(id)
%Excercise 9
%微分方程组 V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t)
%答案 (1)exp(20);(2)0.353;(3)30;(4)451,0.4,9.6
范文三:matlab数学实验答案(胡良剑)
数学实验答案
%Page20,ex1
(5) 等于[exp(1),exp(2);exp(3),exp(4)]
(7) 3=1*3, 8=2*4
(8) a为各列最小值,b为最小值所在的行号
(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture
(11) 答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)
(12) 答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10)
%Page20, ex2
(1)a, b, c的值尽管都是1, 但数据类型分别为数值,字符, 逻辑, 注意a与c相等, 但他们不等于b
(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码
%Page20,ex3
>> r=2;p=0.5;n=12;
>> T=log(r)/n/log(1+0.01*p)
T =
11.5813
%Page20,ex4
>> x=-2:0.05:2;f=x.^4-2.^x;
>> [fmin,min_index]=min(f)
fmin =
-1.3907 %最小值
min_index =
54 %最小值点编址
>> x(min_index)
ans =
0.6500 %最小值点
>> [f1,x1_index]=min(abs(f)) %求近似根--绝对值最小的点
f1 =
0.0328
x1_index =
24
>> x(x1_index)
ans =
-0.8500
>> x(x1_index)=[];f=x.^4-2.^x; %删去绝对值最小的点以求函数绝对值次小的点 >> [f2,x2_index]=min(abs(f)) %求另一近似根--函数绝对值次小的点
f2 =
0.0630
x2_index =
65
>> x(x2_index)
ans =
1.2500
%Page20,ex5
>> z=magic(10)
z =
92 99 1 8 15 67 74 51 58 40
98 80 7 14 16 73 55 57 64 41
4 81 88 20 22 54 56 63 70 47
85 87 19 21 3 60 62 69 71 28
86 93 25 2 9 61 68 75 52 34
17 24 76 83 90 42 49 26 33 65
23 5 82 89 91 48 30 32 39 66
79 6 13 95 97 29 31 38 45 72
10 12 94 96 78 35 37 44 46 53
11 18 100 77 84 36 43 50 27 59
>> sum(z)
ans =
505 505 505 505 505 505 505 505 505 505
>> sum(diag(z))
ans =
505
>> z(:,2)/sqrt(3)
ans =
57.1577
46.1880
46.7654
50.2295
53.6936
13.8564
2.8868
3.4641
6.9282
10.3923
>> z(8,:)=z(8,:)+z(3,:)
z =
92 99 1 8 15 67 74 51 58 40
98 80 7 14 16 73 55 57 64 41
4 81 88 20 22 54 56 63 70 47
85 87 19 21 3 60 62 69 71 28
86 93 25 2 9 61 68 75 52 34
17 24 76 83 90 42 49 26 33 65
23 5 82 89 91 48 30 32 39 66
%Page 40 ex1
先在编辑器窗口写下列M函数,保存为eg2_1.m
function [xbar,s]=ex2_1(x)
n=length(x);
xbar=sum(x)/n;
s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));
例如
>>x=[81 70 65 51 76 66 90 87 61 77];
>>[xbar,s]=ex2_1(x)
xbar =
72.4000
s =
12.1124
%Page 40 ex2
s=log(1);n=0;
while s
n=n+1;
s=s+log(1+n);
end
m=n
计算结果m=37
%Page 40 ex3
clear;
F(1)=1;F(2)=1;k=2;x=0;
e=1e-8; a=(1+sqrt(5))/2;
while abs(x-a)>e
k=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);
end
a,x,k
计算至k=21可满足精度
%Page 40 ex4
clear;tic;s=0;
for i=1:1000000
s=s+sqrt(3)/2^i;
end
s,toc
tic;s=0;i=1;
while i
s=s+sqrt(3)/2^i;i=i+1;
end
s,toc
tic;s=0;
i=1:1000000;
s=sqrt(3)*sum(1./2.^i);
s,toc
%Page 40 ex5
t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16];
plot(t,c)
%Page 40 ex6
%(1)
x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)
y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2])
%(2)参数方法
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t); plot(x,y)
%(3)
x=-3:0.1:3;y=x;
[x,y]=meshgrid(x,y);
z=x.^2+y.^2;
surf(x,y,z)
%(4)
x=-3:0.1:3;y=-3:0.1:13;
[x,y]=meshgrid(x,y);
z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;
surf(x,y,z)
%(5)
t=0:0.01:2*pi;
x=sin(t);y=cos(t);z=cos(2*t);
plot3(x,y,z)
%(6)
theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20);
[theta,fai]=meshgrid(theta,fai);
x=2*sin(fai).*cos(theta);
y=2*sin(fai).*sin(theta);z=2*cos(fai);
surf(x,y,z)
%(7)
x=linspace(0,pi,100);
y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);
plot(x,y1,x,y2,x,y3)
%page41, ex7
x=-1.5:0.05:1.5;
y=1.1*(x>1.1)+x.*(x=-1.1)-1.1*(x
plot(x,y)
%page41,ex8
分别使用which trapz, type trapz, dir C:\MATLAB7\toolbox\matlab\datafun\ %page41,ex9
clear;close;
x=-2:0.1:2;y=x;
[x,y]=meshgrid(x,y);
a=0.5457;b=0.7575;
p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);
p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y
p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y
mesh(x,y,p)
%page41, ex10
lookfor lyapunov
help lyap
>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];
>> X=lyap(A,C)
X =
1.0000 -1.0000 -0.0000
-1.0000 2.0000 1.0000
-0.0000 1.0000 7.0000
Chapter 3
%Exercise 1
>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\b
ans =
0.5000 0.5000 1.0000
ans =
2 2 1
ans =
0.6552 %一元方程组x[2,4,3]=[1,2,3]的近似解
ans =
0 0 0
0 0 0
0.6667 1.3333 1.0000
%矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解
%Exercise 2(1)
>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];
>> rank(A), rank([A,b]) %[A,b]为增广矩阵
ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
2.3830
1.4894
2.0213
%Exercise 2(2)
>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];
>> rank(A), rank([A,b])
ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
-0.4706
-0.2941
%Exercise 2(3)
>> A=[4 1;3 2;1 -5];b=[1;1;1];
>> rank(A), rank([A,b])
ans =
2
ans =
3 %可见方程组无解
>> x=A\b
x =
0.3311
-0.1219 %最小二乘近似解
%Exercise 2(4)
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法 >> rank(a),rank([a,b])
ans =
3
ans =
3 %rank(a)==rank([a,b])
>> a\b
ans =
1
1
0 %一个特解
%Exercise 3
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';
>> x=null(a),x0=a\b
x =
-0.6255
0.6255
-0.2085
0.4170
x0 =
1
1
%通解kx+x0
%Exercise 4
>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95];
>> x1=a*x, x2=a^2*x, x10=a^10*x
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> x0=[0.8 0.2]';
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> [v,e]=eig(a)
v =
0.9806 -0.7071
0.1961 0.7071
e =
1.0000 0
0 0.9400
>> v(:,1)./x
ans =
1.1767
1.1767 %成比例,说明x是最大特征值对应的特征向量
%Exercise 5
%用到公式(3.11)(3.12)
>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]';
>> C=B/diag(x)
C =
0.2400 0.4000 0.0500
0.0900 0.2000 0.0100
0.1200 0.0400 0.0900
>> A=eye(3,3)-C
A =
0.7600 -0.4000 -0.0500
-0.0900 0.8000 -0.0100
-0.1200 -0.0400 0.9100
>> D=[17 17 17]';x=A\D
x =
37.5696
25.7862
24.7690
%Exercise 6(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans =
-94
ans =
0.2553 -0.0213 0.0426
0.1596 -0.1383 -0.2234
0.1809 -0.2234 -0.0532
v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
%Exercise 6(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans =
1
ans =
2.0000 -2.0000 1.0000
1.0000 -1.0000 1.0000
2.0000 -3.0000 2.0000
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774
-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0
0 0 1.0000 - 0.0000i
%Exercise 6(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> det(A),inv(A), [v,d]=eig(A)
ans =
1
ans =
68.0000 -41.0000 -17.0000 10.0000
-41.0000 25.0000 10.0000 -6.0000
-17.0000 10.0000 5.0000 -3.0000
10.0000 -6.0000 -3.0000 2.0000
v =
0.8304 0.0933 0.3963 0.3803
-0.5016 -0.3017 0.6149 0.5286
-0.2086 0.7603 -0.2716 0.5520
0.1237 -0.5676 -0.6254 0.5209
d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
%Exercise 6(4)(以n=5为例)
%关键是矩阵的定义
%方法一(三个for)
n=5;
for i=1:n, a(i,i)=5;end
for i=1:(n-1),a(i,i+1)=6;end
for i=1:(n-1),a(i+1,i)=1;end
a
%方法二(一个for)
n=5;a=zeros(n,n);
a(1,1:2)=[5 6];
for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];end
a(n,[n-1 n])=[1 5];
a
%方法三(不用for)
n=5;a=diag(5*ones(n,1));
b=diag(6*ones(n-1,1));
c=diag(ones(n-1,1));
a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)] %下列计算
>> det(a)
ans =
665
>> inv(a)
ans =
0.3173 -0.5865 1.0286 -1.6241 1.9489
-0.0977 0.4887 -0.8571 1.3534 -1.6241
0.0286 -0.1429 0.5429 -0.8571 1.0286
-0.0075 0.0376 -0.1429 0.4887 -0.5865
0.0015 -0.0075 0.0286 -0.0977 0.3173
>> [v,d]=eig(a)
v =
-0.7843 -0.7843 -0.9237 0.9860 -0.9237
0.5546 -0.5546 -0.3771 -0.0000 0.3771
-0.2614 -0.2614 0.0000 -0.1643 0.0000
0.0924 -0.0924 0.0628 -0.0000 -0.0628
-0.0218 -0.0218 0.0257 0.0274 0.0257
d =
0.7574 0 0 0 0
0 9.2426 0 0 0
0 0 7.4495 0 0
0 0 0 5.0000 0
0 0 0 0 2.5505
%Exercise 7(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a)
v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
>> det(v)
ans =
-0.9255 %v行列式正常, 特征向量线性相关,可对角化
>> inv(v)*a*v %验算
ans =
-3.0527 0.0000 -0.0000
0.0000 3.6760 -0.0000
-0.0000 -0.0000 8.3766
>> [v2,d2]=jordan(a) %也可用jordan
v2 =
0.0798 0.0076 0.9127
0.1886 -0.3141 0.1256
-0.1605 -0.2607 0.4213 %特征向量不同
d2 =
8.3766 0 0
0 -3.0527 - 0.0000i 0
0 0 3.6760 + 0.0000i
>> v2\a*v2
ans =
8.3766 0 0.0000
0.0000 -3.0527 0.0000
0.0000 0.0000 3.6760
>> v(:,1)./v2(:,2) %对应相同特征值的特征向量成比例
ans =
2.4491
2.4491
2.4491
%Exercise 7(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774
-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0
0 0 1.0000 - 0.0000i
>> det(v)
ans =
-5.0566e-028 -5.1918e-017i %v的行列式接近0, 特征向量线性相关,不可对角化 >> [v,d]=jordan(a)
v =
1 0 1
1 0 0
1 -1 0
d =
1 1 0
0 1 1
0 0 1 %jordan标准形不是对角的,所以不可对角化
%Exercise 7(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> [v,d]=eig(A)
v =
0.8304 0.0933 0.3963 0.3803
-0.5016 -0.3017 0.6149 0.5286
-0.2086 0.7603 -0.2716 0.5520
0.1237 -0.5676 -0.6254 0.5209
d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
>> inv(v)*A*v
ans =
0.0102 0.0000 -0.0000 0.0000
0.0000 0.8431 -0.0000 -0.0000
-0.0000 0.0000 3.8581 -0.0000
-0.0000 -0.0000 0 30.2887
%本题用jordan不行, 原因未知
%Exercise 7(4)参考6(4)和7(1), 略
%Exercise 8 只有(3)对称, 且特征值全部大于零, 所以是正定矩阵.
%Exercise 9(1)
>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0]
>> rank(a)
ans =
3
>> rank(a(1:3,:))
ans =
2
>> rank(a([1 2 4],:)) %1,2,4行为最大无关组
ans =
3
>> b=a([1 2 4],:)';c=a([3 5],:)';
>> b\c %线性表示的系数
ans =
0.5000 5.0000
-0.5000 1.0000
0 -5.0000
%Exercise 10
>> a=[1 -2 2;-2 -2 4;2 4 -2]
>> [v,d]=eig(a)
v =
0.3333 0.9339 -0.1293
0.6667 -0.3304 -0.6681
-0.6667 0.1365 -0.7327
d =
-7.0000 0 0
0 2.0000 0
0 0 2.0000
>> v'*v
ans =
1.0000 0.0000 0.0000
0.0000 1.0000 0
0.0000 0 1.0000 %v确实是正交矩阵
%Exercise 11
%设经过6个电阻的电流分别为i1, ..., i6. 列方程组如下
%20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;
%i1=i3+i4;i5=i2+i3;i6=i4+i5;
%计算如下
>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0;
0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0;
0 0 0 0 0 0 -1 -1 1];
>>b=[20 5 0 0 0 0 0 0 0]'; A\b
ans =
13.3453
6.4401
8.5420
3.3274
-1.1807
1.6011
1.7263
0.4204
2.1467
%Exercise 12
>> A=[1 2 3;4 5 6;7 8 0];
>> left=sum(eig(A)), right=sum(trace(A))
left =
6.0000
right =
6
>> left=prod(eig(A)), right=det(A) %原题有错, (-1)^n应删去
left =
27.0000
right =
27
>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3))
fA =
1.0e-012 *
0.0853 0.1421 0.0284
0.1421 0.1421 0
-0.0568 -0.1137 0.1705
>> norm(fA) %f(A)范数接近0
ans =
2.9536e-013
%Exercise 1(1)
roots([1 1 1])
%Exercise 1(2)
roots([3 0 -4 0 2 -1])
%Exercise 1(3)
p=zeros(1,24);
p([1 17 18 22])=[5 -6 8 -5];
roots(p)
%Exercise 1(4)
p1=[2 3];
p2=conv(p1, p1);
p3=conv(p1, p2);
p3(end)=p3(end)-4; %原p3最后一个分量-4
roots(p3)
%Exercise 2
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
fzero(fun,2)
%Exercise 3
fun=inline('x^4-2^x');
fplot(fun,[-2 2]);grid on;
fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)
%Exercise 4
fun=inline('x*sin(1/x)','x');
fplot(fun, [-0.1 0.1]);
x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;
x=[x,-x]
%Exercise 5
fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^2-16*x(3)^2]','x');
[a,b,c]=fsolve(fun,[0 0 0])
%Exercise 6
fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];
[a,b,c]=fsolve(fun,[0.5 0.5])
%Exercise 7
clear; close; t=0:pi/100:2*pi;
x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);
x2=3+sqrt(2)*cos(t); y2=6*sin(t);
plot(x1,y1,x2,y2); grid on; %作图发现4个解的大致位置,然后分别求解
y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])
y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])
y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])
y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])
%Exercise 8(1)
clear;
fun=inline('x.^2.*sin(x.^2-x-2)');
fplot(fun,[-2 2]);grid on; %作图观察
x(1)=-2;
x(3)=fminbnd(fun,-1,-0.5);
x(5)=fminbnd(fun,1,2);
fun2=inline('-x.^2.*sin(x.^2-x-2)');
x(2)=fminbnd(fun2,-2,-1);
x(4)=fminbnd(fun2,-0.5,0.5);
x(6)=2
feval(fun,x)
%答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小, x(2)最大。
%Exercise 8(2)
clear;
fun=inline('3*x.^5-20*x.^3+10');
fplot(fun,[-3 3]);grid on;%作图观察
x(1)=-3;
x(3)=fminsearch(fun,2.5);
fun2=inline('-(3*x.^5-20*x.^3+10)');
x(2)=fminsearch(fun2,-2.5);
x(4)=3;
feval(fun,x)
%Exercise 8(3)
fun=inline('abs(x^3-x^2-x-2)');
fplot(fun,[0 3]);grid on;%作图观察
fminbnd(fun,1.5,2.5)
fun2=inline('-abs(x^3-x^2-x-2)');
fminbnd(fun2,0.5,1.5)
%Exercise 9
close;
x=-2:0.1:1;y=-7:0.1:1;
[x,y]=meshgrid(x,y);
z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9;
mesh(x,y,z);grid on;%作图观察
fun=inline('x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9');
x=fminsearch(fun,[0 0])%求极小值
fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)');
x=fminsearch(fun2,[0 -5])%求极大值
%Exercise 10
clear;t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16];
p2=polyfit(t,c,2)
p3=polyfit(t,c,3)
fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t');
a=lsqcurvefit(fun,[0 0],t,c)%初值可以试探
f=feval(fun, a,t)
norm(f-c)%拟合效果
plot(t,c,t,f) %作图检验
fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改f(x)+20 b=lsqcurvefit(fun2,[0 0],t,c)
figure
f2=feval(fun2, b,t)
norm(f2-c)%拟合效果
plot(t,c,t,f2) %作图检验
%Exercise 11
fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)'); x=fzero(fun, 0, 1)
%Exercise 12
r=5.04/12/100;N=20*12;
x=7500*180 %房屋总价格
y=x*0.3 %首付款额
x0=x-y%贷款总额
a=(1+r)^N*r*x0/((1+r)^N-1)%月付还款额
r1=4.05/12/100;x1=10*10000;%公积金贷款
a1=(1+r1)^N*r1*x1/((1+r1)^N-1)
x2=x0-x1%商业贷款
a2=(1+r)^N*r*x2/((1+r)^N-1)
a=a1+a2
%Exercise 13
%列方程th*R^2+(pi-2*th)*r^2-R*r*sin(th)=pi*r^2/2
%化简得sin(2*th)-2*th*cos(2*th)=pi/2
%以下Matlab计算
clear;fun= inline('sin(2*th)-2*th*cos(2*th)-pi/2','th')
th=fsolve(fun,pi/4)
R=20*cos(th)
%Exercise 14
%先在Editor窗口写M函数保存
function x=secant(fname,x0,x1,e)
while abs(x0-x1)>e,
x=x1-(x1-x0)*feval(fname,x1)/(feval(fname,x1)-feval(fname,x0));
x0=x1;x1=x;
end
%再在指令窗口
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
secant(fun,1,2,1e-8)
%Exercise 15
%作系数为a,初值为xo,从第m步到第n步迭代过程的M函数: function f=ex4_15fun(a,x0,m,n)
x(1)=x0; y(1)=a*x(1)+1;x(2)=y(1);
if m
for i=2:n
y(i)=a*x(i)+1; x(i+1)=y(i);
if i>m, plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]); end
end
hold off;
%M脚本文件
subplot(2,2,1);ex4_15fun(0.9,1,1,20);
subplot(2,2,2);ex4_15fun(-0.9,1,1,20);
subplot(2,2,3);ex4_15fun(1.1,1,1,20);
subplot(2,2,4);ex4_15fun(-1.1,1,1,20);
%Exercise 16
%设夹角t, 问题转化为 min f=5/sin(t)+10/cos(t)
%取初始值pi/4, 计算如下
fun=@(t)5/sin(t)+10/cos(t);
[t,f]=fminsearch(fun, pi/4)
t =
0.6709
f =
20.8097
%Exercise 17
%提示:x(k+2)=f(x(k))=a^2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k))) %计算平衡点x
%|f'(x)|
%Exercise 18
%先写M文件
function f=ex4_18(a,x0,n)
x=zeros(1,n);y=x;
x(1)=x0;
y(1)=a*x(1)+1;
x(2)=y(1);
plot([x(1),x(1),x(2)],[0,y(1),y(1)],'r');
hold on;
for i=2:n
y(i)=a*x(i)+1;
x(i+1)=y(i);
plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]) end
hold off;
%再执行指令
>> ex4_18(0.9,1,20)
>> ex4_18(-0.9,1,20)
>> ex4_18(1.1,1,20)
>> ex4_18(-1.1,1,20)
%Exercise 19
clear; close; x(1)=0; y(1)=0;
for k=1:3000
x(k+1)=1+y(k)-1.4*x(k)^2; y(k+1)=0.3*x(k); end
plot(x(1000:1500),y(1000:1500),'+g');hold on plot(x(1501:2000),y(1501:2000),'.b'); plot(x(2001:2500),y(2001:2500),'*y'); plot(x(2501:3001),y(2501:3001),'.r');
%Exercise 1
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0];
trapz(x,y)
%Exercise 2
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0];
diff(y)./diff(x)
%Exercise 3
xa=-1:0.1:1;ya=0:0.1:2;
[x,y]=meshgrid(xa,ya);
z=x.*exp(-x.^2 -y.^3);
[px,py] = gradient(z,xa,ya);
px
%Exercise 4
t=0:0.01:1.5;
x=log(cos(t));
y=cos(t)-t.*sin(t);
dydx=gradient(y,x)
[x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点
dydx(id)
%Exercise 5(2)
fun=inline('exp(2*x).*cos(x).^3');
quadl(fun,0,2*pi)
或用trapz
x=linspace(0,2*pi,100);
y=exp(2*x).*cos(x).^3;
trapz(x,y)
%Exercise 5(3)
fun=@(x)x.*log(x.^4).*asin(1./x.^2);
quadl(fun,1,3)
或用trapz
x=1:0.01:3;
y=feval(fun,x);
trapz(x,y)
%Exercise 5(4)
fun=@(x)sin(x)./x;
quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替 %Exercise 5(5)
%参考Exercise 5(4)
%Exercise 5(6)
fun=inline('sqrt(1+r.^2.*sin(th))','r','th');
dblquad(fun,0,1,0,2*pi)
%Exercise 5(7)
首先建立84页函数dblquad2
clear;
fun=@(x,y)1+x+y.^2;
clo=@(x)-sqrt(2*x-x.^2);
dup=@(x)sqrt(2*x-x.^2);
dblquad2(fun,0,2,clo,dhi,100)
%Exercise 6
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t);
dx=gradient(x,t);dy=gradient(y,t);
f=sqrt(dx.^2+dy.^2);
trapz(t,f)
%Exercise 7
xa=-1:0.1:1;ya=0:0.1:2;
[x,y]=meshgrid(xa,ya);
z=x.*exp(x.^2+y.^2);
[zx,zy]=gradient(z,xa,ya);
f=sqrt(1+zx.^2+zy.^2);
s=0;
for i=2:length(xa)
for j=2:length(ya)
s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4;
end
end
s
%Exercise 8
funl=inline('-(-x).^0.2.*cos(x)');
funr=inline('x.^0.2.*cos(x)');
quadl(funl,-1,0)+quadl(funr,0,1)
%Exercise 9 (以I32为例)
fun=@(x)abs(sin(x));
h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y)
h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致,结果失真
q1=quad(fun,0,32*pi)
q2=quadl(fun,0,32*pi)
%Exercise 10(2)
%先在程序编辑器,写下列函数,保存为ex5_10_2f
function d=ex5_10_2f(fname,a,h0,e)
h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);
d0=d+2*e;
while abs(d-d0)>e
d0=d;h0=h;h=h0/2;
d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);
end
%再在指令窗口执行
fun=inline('x.^2*sin(x.^2-x-2)','x');
d=ex5_10_2f(fun,1.4,0.1,1e-3)
%Exercise 11
%提示:f上升时,f'>0;f下降时,f'
%Exercise 12
在程序编辑器,写下列函数,保存为ex5_12f
function I=ex5_12(fname,a,b,n)
x=linspace(a,b,n+1);
y=feval(fname,x);
I=(b-a)/n/3*(y(1)+y(n+1)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));
%再在指令窗口执行
ex5_12(inline('1/sqrt(2*pi)*exp(-x.^2/2)'),0,1,50)
%Exercise 13
fun=inline('5400*v./(8.276*v.^2+2000)','v');
quadl(fun,15,30)
%Exercise 14
重心不超过凳边沿。1/2, 2/3, 3/4, ...,n/(n+1)
%Exercise15
利润函数fun=inline('(p-c0+k*log(M*exp(-a*p)))*M*exp(-a*p)','p');
求p使fun最大
%Exercise 16
clear; x=-3/4:0.01:3/4;
y=(3/4+x)*2.*sqrt(1-16/9.*x.^2)*9.8;
P=trapz(x,y) %单位:千牛
%Exercise 17
clear; close;
fplot('17-t^(2/3)-5-2*t^(2/3)',[0,20]); grid;
t=fzero('17-x^(2/3)-5-2*x^(2/3)',7)
t=0:0.1:8; y=17-t.^(2/3)-5-2*t.^(2/3);
trapz(t,y)-20 %单位:百万元
%Exercise 18
%曲面面积计算
%Excercise 1(1)
fun=inline('x+y','x','y');
[t,y]=ode45(fun,[0 1 2 3],1) %注意由于初值为y(0)=1,[0 1 2 3]中0不可缺
%Excercise 1(3)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=0.01*y(2)^2-2*y(1)+sin(t)
%运行下列指令
clear;close;
fun=@(t,y)[y(2);0.01*y(2)^2-2*y(1)+sin(t)];
[t,y]=ode45(fun,[0 5],[0;1]);
plot(t,y(:,1))
%Excercise 1(5)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=-mu*(y(1)^2-1)*y(2)-y(1)
%运行下列指令,注意参数mu的处理
clear;close;
fun=@(t,y,mu)[y(2);-mu*(y(1)^2-1)*y(2)-y(1)];
[t,y]=ode45(fun,[0 20],[2;0],[],1);
plot(y(:,1),y(:,2));hold on;
[t,y]=ode45(fun,[0 20],[2;0],[],2);
plot(y(:,1),y(:,2),'r');hold off;
%Excercise 2
roots([1 10 54 132 137 50])
%通解A1*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-t)+A5*t*exp(-t) %Excercise 3
dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');
[x,y]=ode45(dfun,[0,50],[1;-1]);length(x)
%所用节点很多
[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x)
%所用节点很少
%Excercise 4
clear;
dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3;
x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3]','t','x');
[t,x]=ode45(dfun,[0 24],[1.2; 0; 0; -1.04935371]);
plot(x(:,1),x(:,3));
%Excercise 5
%方程y'=2x+y^2,y(0)=0
clear;close;
fun=inline('2*x+y^2','x','y');
[x,y]=ode45(fun,[0 1.57],0); %x的上界再增加,解会
plot(x,y)
%Excercise 6
clear;close;
fun=@(t,x,a,b)a*x+b;
[t,x]=ode45(fun,[0 10],0.1,[],1,1);
subplot(2,4,1);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,1);
subplot(2,4,2);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],1,-1);
subplot(2,4,3);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,-1);
subplot(2,4,4);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,1);
subplot(2,4,5);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,1);
subplot(2,4,6);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,-1);
subplot(2,4,7);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,-1);
subplot(2,4,8);plot(t,x)
%Excercise 7
%微分方程 T'=k(c-T),T(0)=20
dsolve('DT=k*(c-T)','T(0)=20','t')
%得c+exp(-k*t)*(-c+20)
%利用T(10)=25.2, T(20)=28.32拟合(或者解非线性方程)
fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t')
lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32])
%解得户外温度c=33,比例系数k=0.05.
%Excercise 8
%微分方程 x'/x=0.5*(1-x),x(0)=0.05
fun=inline('0.5*(1-x)*x','t','x');
[t,x]=ode45(fun,[0 10],0.05);
plot(t,x)
id=min(find(x>0.8));
t(id)
%Excercise 9
%微分方程组 V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t)
%答案(1)exp(20);(2)0.353;(3)30;(4)451,0.4,9.6
%Chapter 7
%Exercise 1
syms ph th;
a=sin(ph)*cos(th)-cos(ph)*sin(th)-sin(ph-th);
simple(a)
%化简后差的结果为0
%Exercise 2
syms x;s=x^4-5*x^3+5*x^2+5*x-6;
factor(s)
%Exercise 3
syms a;A=[1 2;2 a];
iA=inv(A),[v,d]=eig(A)
%Exercise 4
syms x y;
limit((3^x+9^x)^(1/x),x,inf)
s1=limit(x*y/(sqrt(x*y+1)-1),x,0);s2=limit(s1,y,0)
%Exercise 5
syms k n x;s1=symsum(k^2,k,1,n);s1=simple(s1)
s2=symsum(k^(-2),k,1,inf);s2=simple(s2)
s3=symsum(1/(2*n+1)/(2*x+1)^(2*n+1),n,0,inf);s3=simple(s3) %Exercise 6
syms x y z;s=sin(x^2*y*z);
s=diff(s,x,2);
s=diff(s,y,1);
s=subs(s,{x,y,z},{1,1,3})
%Exercise 7
syms x;s=log(x+sqrt(1+x^2));taylor(s,8,0,x)
%Exercise 8 (以第四章习题9为例)
先用符号运算求偏导数
syms x y;f=y^3/9+3*x^2*y+9*x^2+y^2+x*y+9;
fx=diff(f,x),fy=diff(f,y)
根据计算结果得方程组. 求解方程组
[sx,sy]=solve(6*x*y+18*x+y,1/3*y^2+3*x^2+2*y+x,x,y)
得四个解(0,0),(-1/3,-6),(-7/6,-7/2),(5/6,-5/2).计算Hesse矩阵 fh2=[diff(fx,x),diff(fx,y);diff(fy,x),diff(fy,y)]
计算
eig(subs(fh2,[x,y],[0,0]))
得知正定,所以是极小值点.极小值用
subs(f,[x,y],[0,0])
求得。同理可得(-1/3,-6)为极大值点,其它两个为鞍点。 %Exercise 9(以第一小题为例)
syms y;f=exp(2*y)/(exp(y)+2);
fi=int(f,y)
s=simple(diff(fi)-f)
%Exercise 10
syms x y;f=(x-y)^3*sin(x+2*y);Ix=simple(int(f,y,-x,x))
%Exercise 11(3)
syms x;f=x*log(x^4)*asin(1/x);Ix=int(f,x,1,3);vpa(Ix)
%Exercise 12
%1(3)
syms x;solve(5*x^23-6*x^7+8*x^6-5*x^2)
%6
syms a b;s=solve(a-0.7*sin(a)-0.2*cos(b),b-0.7*cos(a)+0.2*sin(b)); s.a,s.b
%Exercise 13
%1(3)
dsolve('D2y-0.01*Dy^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')(解不出) %1(4)
dsolve('2*D2x-5*Dx+3*x=45*exp(2*t)','x(0)=2','Dx(0)=1','t') %Exercise 14
%6(ii)
ezplot('x^2/4+y^2/9=1')
%6(vi)
ezmesh('2*sin(ph)*cos(th)','2*sin(ph)*sin(th)','2*cos(ph)',[0 pi/2 0 2*pi])
范文四:matlab数学实验答案(胡良剑)版本
%Page20,ex1
(5) 等于[exp(1),exp(2);exp(3),exp(4)]
(7) 3=1*3, 8=2*4
(8) a为各列最小值,b为最小值所在的行号
(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture
(11) 答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)
(12) 答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10)
%Page20, ex2
(1)a, b, c的值尽管都是1, 但数据类型分别为数值,字符, 逻辑, 注意a与c相等,他们不等于b
(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码
%Page20,ex3
>> r=2;p=0.5;n=12;
>> T=log(r)/n/log(1+0.01*p)
T =
11.5813
%Page20,ex4
>> x=-2:0.05:2;f=x.^4-2.^x;
>> [fmin,min_index]=min(f)
fmin =
-1.3907 %最小值
min_index =
54 %最小值点编址
>> x(min_index)
ans =
0.6500 %最小值点
>> [f1,x1_index]=min(abs(f)) %求近似根--绝对值最小的点
f1 =
0.0328
x1_index =
24
>> x(x1_index)
ans =
-0.8500
>> x(x1_index)=[];f=x.^4-2.^x; %删去绝对值最小的点以求函数绝对值次小的点 >> [f2,x2_index]=min(abs(f)) %求另一近似根--函数绝对值次小的点
f2 =
0.0630
x2_index =
65 但
>> x(x2_index)
ans =
1.2500
%Page20,ex5
>> z=magic(10)
z =
92 99 1 8 15 67 74 51 58 40
98 80 7 14 16 73 55 57 64 41
4 81 88 20 22 54 56 63 70 47
85 87 19 21 3 60 62 69 71 28
86 93 25 2 9 61 68 75 52 34
17 24 76 83 90 42 49 26 33 65
23 5 82 89 91 48 30 32 39 66
79 6 13 95 97 29 31 38 45 72
10 12 94 96 78 35 37 44 46 53
11 18 100 77 84 36 43 50 27 59
>> sum(z)
ans =
505 505 505 505 505 505 505 505 505 505
>> sum(diag(z))
ans =
505
>> z(:,2)/sqrt(3)
ans =
57.1577
46.1880
46.7654
50.2295
53.6936
13.8564
2.8868
3.4641
6.9282
10.3923
>> z(8,:)=z(8,:)+z(3,:)
z =
92 99 1 8 15 67 74 51 58 40
98 80 7 14 16 73 55 57 64 41
4 81 88 20 22 54 56 63 70 47
85 87 19 21 3 60 62 69 71 28
86 93 25 2 9 61 68 75 52 34
17 24 76 83 90 42 49 26 33 65
23 5 82 89 91 48 30 32 39 66
83 87 101 115 119 83 87 101 115 119
10 12 94 96 78 35 37 44 46 53
11 18 100 77 84 36 43 50 27 59
%Page 40 ex1
先在编辑器窗口写下列M函数,保存为eg2_1.m
function [xbar,s]=ex2_1(x)
n=length(x);
xbar=sum(x)/n;
s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));
例如
>>x=[81 70 65 51 76 66 90 87 61 77];
>>[xbar,s]=ex2_1(x)
xbar =
72.4000
s =
12.1124
%Page 40 ex2
s=log(1);n=0;
while s
n=n+1;
s=s+log(1+n);
end
m=n
计算结果m=37
%Page 40 ex3
clear;
F(1)=1;F(2)=1;k=2;x=0;
e=1e-8; a=(1+sqrt(5))/2;
while abs(x-a)>e
k=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);
end
a,x,k
计算至k=21可满足精度
%Page 40 ex4
clear;tic;s=0;
for i=1:1000000
s=s+sqrt(3)/2^i;
end
s,toc
tic;s=0;i=1;
while i
s=s+sqrt(3)/2^i;i=i+1;
end
s,toc
tic;s=0;
i=1:1000000;
s=sqrt(3)*sum(1./2.^i);
s,toc
%Page 40 ex5
t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16];
plot(t,c)
%Page 40 ex6
%(1)
x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)
y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2])
%(2)参数方法
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t); plot(x,y)
%(3)
x=-3:0.1:3;y=x;
[x,y]=meshgrid(x,y);
z=x.^2+y.^2;
surf(x,y,z)
%(4)
x=-3:0.1:3;y=-3:0.1:13;
[x,y]=meshgrid(x,y);
z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;
surf(x,y,z)
%(5)
t=0:0.01:2*pi;
x=sin(t);y=cos(t);z=cos(2*t);
plot3(x,y,z)
%(6)
theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20);
[theta,fai]=meshgrid(theta,fai);
x=2*sin(fai).*cos(theta);
y=2*sin(fai).*sin(theta);z=2*cos(fai);
surf(x,y,z)
%(7)
x=linspace(0,pi,100);
y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);
plot(x,y1,x,y2,x,y3)
%page41, ex7
x=-1.5:0.05:1.5;
y=1.1*(x>1.1)+x.*(x=-1.1)-1.1*(x
plot(x,y)
%page41,ex8
分别使用which trapz, type trapz, dir C:\MATLAB7\toolbox\matlab\datafun\ %page41,ex9
clear;close;
x=-2:0.1:2;y=x;
[x,y]=meshgrid(x,y);
a=0.5457;b=0.7575;
p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);
p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y
p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y
mesh(x,y,p)
%page41, ex10
lookfor lyapunov
help lyap
>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];
>> X=lyap(A,C)
X =
1.0000 -1.0000 -0.0000
-1.0000 2.0000 1.0000
-0.0000 1.0000 7.0000
%Chapter 3
%Exercise 1
>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\b
ans =
0.5000 0.5000 1.0000
ans =
2 2 1
ans =
0.6552 %一元方程组x[2,4,3]=[1,2,3]的近似解
ans =
0 0 0
0 0 0
0.6667 1.3333 1.0000
%矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解 %Exercise 2(1)
>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];
>> rank(A), rank([A,b]) %[A,b]为增广矩阵
ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
2.3830
1.4894
2.0213
%Exercise 2(2)
>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];
>> rank(A), rank([A,b])
ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
-0.4706
-0.2941
%Exercise 2(3)
>> A=[4 1;3 2;1 -5];b=[1;1;1];
>> rank(A), rank([A,b])
ans =
2
ans =
3 %可见方程组无解
>> x=A\b
x =
0.3311
-0.1219 %最小二乘近似解
%Exercise 2(4)
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法 >> rank(a),rank([a,b])
ans =
3
ans =
3 %rank(a)==rank([a,b])
>> a\b
ans =
1
1
0 %一个特解
%Exercise 3
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';
>> x=null(a),x0=a\b
x =
-0.6255
0.6255
-0.2085
0.4170
x0 =
1
1
%通解kx+x0
%Exercise 4
>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95];
>> x1=a*x, x2=a^2*x, x10=a^10*x
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> x0=[0.8 0.2]';
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> [v,e]=eig(a)
v =
0.9806 -0.7071
0.1961 0.7071
e =
1.0000 0
0 0.9400
>> v(:,1)./x
ans =
1.1767
1.1767 %成比例,说明x是最大特征值对应的特征向量 %Exercise 5
%用到公式(3.11)(3.12)
>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]';
>> C=B/diag(x)
C =
0.2400 0.4000 0.0500
0.1200 0.0400 0.0900
>> A=eye(3,3)-C
A =
0.7600 -0.4000 -0.0500
-0.0900 0.8000 -0.0100
-0.1200 -0.0400 0.9100
>> D=[17 17 17]';x=A\D
x =
37.5696
25.7862
24.7690
%Exercise 6(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a) ans =
-94
ans =
0.2553 -0.0213 0.0426
0.1596 -0.1383 -0.2234
0.1809 -0.2234 -0.0532
v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
%Exercise 6(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a) ans =
1
ans =
2.0000 -2.0000 1.0000
1.0000 -1.0000 1.0000
2.0000 -3.0000 2.0000
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774
-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0
%Exercise 6(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> det(A),inv(A), [v,d]=eig(A)
ans =
1
ans =
68.0000 -41.0000 -17.0000 10.0000
-41.0000 25.0000 10.0000 -6.0000
-17.0000 10.0000 5.0000 -3.0000
10.0000 -6.0000 -3.0000 2.0000
v =
0.8304 0.0933 0.3963 0.3803
-0.5016 -0.3017 0.6149 0.5286
-0.2086 0.7603 -0.2716 0.5520
0.1237 -0.5676 -0.6254 0.5209
d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
%Exercise 6(4)(以n=5为例)
%关键是矩阵的定义
%方法一(三个for)
n=5;
for i=1:n, a(i,i)=5;end
for i=1:(n-1),a(i,i+1)=6;end
for i=1:(n-1),a(i+1,i)=1;end
a
%方法二(一个for)
n=5;a=zeros(n,n);
a(1,1:2)=[5 6];
for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];end
a(n,[n-1 n])=[1 5];
a
%方法三(不用for)
n=5;a=diag(5*ones(n,1));
b=diag(6*ones(n-1,1));
c=diag(ones(n-1,1));
a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)] %下列计算
>> det(a)
ans =
665
>> inv(a)
ans =
0.3173 -0.5865 1.0286 -1.6241 1.9489
-0.0977 0.4887 -0.8571 1.3534 -1.6241
0.0286 -0.1429 0.5429 -0.8571 1.0286
-0.0075 0.0376 -0.1429 0.4887 -0.5865
0.0015 -0.0075 0.0286 -0.0977 0.3173
>> [v,d]=eig(a)
v =
-0.7843 -0.7843 -0.9237 0.9860 -0.9237
0.5546 -0.5546 -0.3771 -0.0000 0.3771
-0.2614 -0.2614 0.0000 -0.1643 0.0000
0.0924 -0.0924 0.0628 -0.0000 -0.0628
-0.0218 -0.0218 0.0257 0.0274 0.0257
d =
0.7574 0 0 0 0
0 9.2426 0 0 0
0 0 7.4495 0 0
0 0 0 5.0000 0
0 0 0 0 2.5505
%Exercise 7(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a)
v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
>> det(v)
ans =
-0.9255 %v行列式正常, 特征向量线性相关,可对角化 >> inv(v)*a*v %验算
ans =
-3.0527 0.0000 -0.0000
0.0000 3.6760 -0.0000
-0.0000 -0.0000 8.3766
>> [v2,d2]=jordan(a) %也可用jordan
v2 =
0.0798 0.0076 0.9127
0.1886 -0.3141 0.1256
-0.1605 -0.2607 0.4213 %特征向量不同
d2 =
8.3766 0 0
0 -3.0527 - 0.0000i 0
0 0 3.6760 + 0.0000i
>> v2\a*v2
ans =
8.3766 0 0.0000
0.0000 -3.0527 0.0000
0.0000 0.0000 3.6760
>> v(:,1)./v2(:,2) %对应相同特征值的特征向量成比例
ans =
2.4491
2.4491
2.4491
%Exercise 7(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774
-0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0
0 0 1.0000 - 0.0000i
>> det(v)
ans =
-5.0566e-028 -5.1918e-017i %v的行列式接近0, 特征向量线性相关,不可对角化 >> [v,d]=jordan(a)
v =
1 0 1
1 0 0
1 -1 0
d =
1 1 0
0 1 1
0 0 1 %jordan标准形不是对角的,所以不可对角化
%Exercise 7(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> [v,d]=eig(A)
v =
0.8304 0.0933 0.3963 0.3803
-0.5016 -0.3017 0.6149 0.5286
-0.2086 0.7603 -0.2716 0.5520
0.1237 -0.5676 -0.6254 0.5209
d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
>> inv(v)*A*v
ans =
0.0102 0.0000 -0.0000 0.0000
0.0000 0.8431 -0.0000 -0.0000
-0.0000 0.0000 3.8581 -0.0000
-0.0000 -0.0000 0 30.2887
%本题用jordan不行, 原因未知
%Exercise 7(4)参考6(4)和7(1), 略
%Exercise 8 只有(3)对称, 且特征值全部大于零, 所以是正定矩阵.
%Exercise 9(1)
>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0]
>> rank(a)
ans =
3
>> rank(a(1:3,:))
ans =
2
>> rank(a([1 2 4],:)) %1,2,4行为最大无关组
ans =
3
>> b=a([1 2 4],:)';c=a([3 5],:)';
>> b\c %线性表示的系数
ans =
0.5000 5.0000
-0.5000 1.0000
0 -5.0000
%Exercise 10
>> a=[1 -2 2;-2 -2 4;2 4 -2]
>> [v,d]=eig(a)
v =
0.3333 0.9339 -0.1293
0.6667 -0.3304 -0.6681
-0.6667 0.1365 -0.7327
d =
-7.0000 0 0
0 2.0000 0
0 0 2.0000
>> v'*v
ans =
1.0000 0.0000 0.0000
0.0000 1.0000 0
0.0000 0 1.0000 %v确实是正交矩阵
%Exercise 11
%设经过6个电阻的电流分别为i1, ..., i6. 列方程组如下
%20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;
%i1=i3+i4;i5=i2+i3;i6=i4+i5;
%计算如下
>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0;
0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0;
0 0 0 0 0 0 -1 -1 1];
>>b=[20 5 0 0 0 0 0 0 0]'; A\b
ans =
13.3453
6.4401
8.5420
3.3274
-1.1807
1.6011
1.7263
0.4204
2.1467
%Exercise 12
>> A=[1 2 3;4 5 6;7 8 0];
>> left=sum(eig(A)), right=sum(trace(A))
left =
6.0000
right =
6
>> left=prod(eig(A)), right=det(A) %原题有错, (-1)^n应删去
left =
27.0000
right =
27
>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3))
fA =
1.0e-012 *
0.0853 0.1421 0.0284
0.1421 0.1421 0
-0.0568 -0.1137 0.1705
>> norm(fA) %f(A)范数接近0
ans =
2.9536e-013
%Exercise 1(1)
roots([1 1 1])
%Exercise 1(2)
roots([3 0 -4 0 2 -1])
%Exercise 1(3)
p=zeros(1,24);
p([1 17 18 22])=[5 -6 8 -5];
roots(p)
%Exercise 1(4)
p1=[2 3];
p2=conv(p1, p1);
p3=conv(p1, p2);
p3(end)=p3(end)-4; %原p3最后一个分量-4
roots(p3)
%Exercise 2
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
fzero(fun,2)
%Exercise 3
fun=inline('x^4-2^x');
fplot(fun,[-2 2]);grid on;
fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)
%Exercise 4
fun=inline('x*sin(1/x)','x');
fplot(fun, [-0.1 0.1]);
x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;
x=[x,-x]
%Exercise 5
fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^2-16*x(3)^2]','x');
[a,b,c]=fsolve(fun,[0 0 0])
%Exercise 6
fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];
[a,b,c]=fsolve(fun,[0.5 0.5])
%Exercise 7
clear; close; t=0:pi/100:2*pi;
x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);
x2=3+sqrt(2)*cos(t); y2=6*sin(t);
plot(x1,y1,x2,y2); grid on; %作图发现4个解的大致位置,然后分别求解
y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])
y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])
y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])
y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])
%Exercise 8(1)
clear;
fun=inline('x.^2.*sin(x.^2-x-2)');
fplot(fun,[-2 2]);grid on; %作图观察
x(1)=-2;
x(3)=fminbnd(fun,-1,-0.5);
x(5)=fminbnd(fun,1,2);
fun2=inline('-x.^2.*sin(x.^2-x-2)');
x(2)=fminbnd(fun2,-2,-1);
x(4)=fminbnd(fun2,-0.5,0.5);
x(6)=2
feval(fun,x)
%答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小, x(2)最大。
%Exercise 8(2)
clear;
fun=inline('3*x.^5-20*x.^3+10');
fplot(fun,[-3 3]);grid on;%作图观察
x(1)=-3;
x(3)=fminsearch(fun,2.5);
fun2=inline('-(3*x.^5-20*x.^3+10)');
x(2)=fminsearch(fun2,-2.5);
x(4)=3;
feval(fun,x)
%Exercise 8(3)
fun=inline('abs(x^3-x^2-x-2)');
fplot(fun,[0 3]);grid on;%作图观察
fminbnd(fun,1.5,2.5)
fun2=inline('-abs(x^3-x^2-x-2)');
fminbnd(fun2,0.5,1.5)
%Exercise 9
close;
x=-2:0.1:1;y=-7:0.1:1;
[x,y]=meshgrid(x,y);
z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9;
mesh(x,y,z);grid on;%作图观察
fun=inline('x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9');
x=fminsearch(fun,[0 0])%求极小值
fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)');
x=fminsearch(fun2,[0 -5])%求极大值
%Exercise 10
clear;t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16];
p2=polyfit(t,c,2)
p3=polyfit(t,c,3)
fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t');
a=lsqcurvefit(fun,[0 0],t,c)%初值可以试探
f=feval(fun, a,t)
norm(f-c)%拟合效果
plot(t,c,t,f) %作图检验
fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改f(x)+20
b=lsqcurvefit(fun2,[0 0],t,c)
figure
f2=feval(fun2, b,t)
norm(f2-c)%拟合效果
plot(t,c,t,f2) %作图检验
%Exercise 11
fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)');
x=fzero(fun, 0, 1)
%Exercise 12
r=5.04/12/100;N=20*12;
x=7500*180 %房屋总价格
y=x*0.3 %首付款额
x0=x-y%贷款总额
a=(1+r)^N*r*x0/((1+r)^N-1)%月付还款额
r1=4.05/12/100;x1=10*10000;%公积金贷款
a1=(1+r1)^N*r1*x1/((1+r1)^N-1)
x2=x0-x1%商业贷款
a2=(1+r)^N*r*x2/((1+r)^N-1)
a=a1+a2
%Exercise 13
%列方程th*R^2+(pi-2*th)*r^2-R*r*sin(th)=pi*r^2/2
%化简得sin(2*th)-2*th*cos(2*th)=pi/2
%以下Matlab计算
clear;fun= inline('sin(2*th)-2*th*cos(2*th)-pi/2','th')
th=fsolve(fun,pi/4)
R=20*cos(th)
%Exercise 14
%先在Editor窗口写M函数保存
function x=secant(fname,x0,x1,e)
while abs(x0-x1)>e,
x=x1-(x1-x0)*feval(fname,x1)/(feval(fname,x1)-feval(fname,x0));
x0=x1;x1=x;
end
%再在指令窗口
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
secant(fun,1,2,1e-8)
%Exercise 15
%作系数为a,初值为xo,从第m步到第n步迭代过程的M函数:
function f=ex4_15fun(a,x0,m,n)
x(1)=x0; y(1)=a*x(1)+1;x(2)=y(1);
if m
for i=2:n
y(i)=a*x(i)+1; x(i+1)=y(i);
if i>m, plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]); end
end
hold off;
%M脚本文件
subplot(2,2,1);ex4_15fun(0.9,1,1,20);
subplot(2,2,2);ex4_15fun(-0.9,1,1,20);
subplot(2,2,3);ex4_15fun(1.1,1,1,20);
subplot(2,2,4);ex4_15fun(-1.1,1,1,20);
%Exercise 16
%设夹角t, 问题转化为 min f=5/sin(t)+10/cos(t)
%取初始值pi/4, 计算如下
fun=@(t)5/sin(t)+10/cos(t);
[t,f]=fminsearch(fun, pi/4)
t =
0.6709
f =
20.8097
%Exercise 17
%提示:x(k+2)=f(x(k))=a^2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k)))
%计算平衡点x
%|f'(x)|
%Exercise 18
%先写M文件
function f=ex4_18(a,x0,n)
x=zeros(1,n);y=x;
x(1)=x0;
y(1)=a*x(1)+1;
x(2)=y(1);
plot([x(1),x(1),x(2)],[0,y(1),y(1)],'r');
hold on;
for i=2:n
y(i)=a*x(i)+1;
x(i+1)=y(i);
plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)])
end
hold off;
%再执行指令
>> ex4_18(0.9,1,20)
>> ex4_18(-0.9,1,20)
>> ex4_18(1.1,1,20)
>> ex4_18(-1.1,1,20)
%Exercise 19
clear; close; x(1)=0; y(1)=0;
for k=1:3000
x(k+1)=1+y(k)-1.4*x(k)^2; y(k+1)=0.3*x(k);
end
plot(x(1000:1500),y(1000:1500),'+g');hold on
plot(x(1501:2000),y(1501:2000),'.b');
plot(x(2001:2500),y(2001:2500),'*y');
plot(x(2501:3001),y(2501:3001),'.r');
%Exercise 1
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0];
trapz(x,y)
%Exercise 2
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0];
diff(y)./diff(x)
%Exercise 3
xa=-1:0.1:1;ya=0:0.1:2;
[x,y]=meshgrid(xa,ya);
z=x.*exp(-x.^2 -y.^3);
[px,py] = gradient(z,xa,ya);
px
%Exercise 4
t=0:0.01:1.5;
x=log(cos(t));
y=cos(t)-t.*sin(t);
dydx=gradient(y,x)
[x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点
dydx(id)
%Exercise 5(2)
fun=inline('exp(2*x).*cos(x).^3');
quadl(fun,0,2*pi)
或用trapz
x=linspace(0,2*pi,100);
y=exp(2*x).*cos(x).^3;
trapz(x,y)
%Exercise 5(3)
fun=@(x)x.*log(x.^4).*asin(1./x.^2);
quadl(fun,1,3)
或用trapz
x=1:0.01:3;
y=feval(fun,x);
trapz(x,y)
%Exercise 5(4)
fun=@(x)sin(x)./x;
quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替 %Exercise 5(5)
%参考Exercise 5(4)
%Exercise 5(6)
fun=inline('sqrt(1+r.^2.*sin(th))','r','th');
dblquad(fun,0,1,0,2*pi)
%Exercise 5(7)
首先建立84页函数dblquad2
clear;
fun=@(x,y)1+x+y.^2;
clo=@(x)-sqrt(2*x-x.^2);
dup=@(x)sqrt(2*x-x.^2);
dblquad2(fun,0,2,clo,dhi,100)
%Exercise 6
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t);
dx=gradient(x,t);dy=gradient(y,t);
f=sqrt(dx.^2+dy.^2);
trapz(t,f)
%Exercise 7
xa=-1:0.1:1;ya=0:0.1:2;
[x,y]=meshgrid(xa,ya);
z=x.*exp(x.^2+y.^2);
[zx,zy]=gradient(z,xa,ya);
f=sqrt(1+zx.^2+zy.^2);
s=0;
for i=2:length(xa)
for j=2:length(ya)
s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4; end
end
s
%Exercise 8
funl=inline('-(-x).^0.2.*cos(x)');
funr=inline('x.^0.2.*cos(x)');
quadl(funl,-1,0)+quadl(funr,0,1)
%Exercise 9 (以I32为例)
fun=@(x)abs(sin(x));
h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y)
h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致,结果失真 q1=quad(fun,0,32*pi)
q2=quadl(fun,0,32*pi)
%Exercise 10(2)
%先在程序编辑器,写下列函数,保存为ex5_10_2f
function d=ex5_10_2f(fname,a,h0,e)
h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h); d0=d+2*e;
while abs(d-d0)>e
d0=d;h0=h;h=h0/2;
d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);
end
%再在指令窗口执行
fun=inline('x.^2*sin(x.^2-x-2)','x');
d=ex5_10_2f(fun,1.4,0.1,1e-3)
%Exercise 11
%提示:f上升时,f'>0;f下降时,f'
%Exercise 12
在程序编辑器,写下列函数,保存为ex5_12f
function I=ex5_12(fname,a,b,n)
x=linspace(a,b,n+1);
y=feval(fname,x);
I=(b-a)/n/3*(y(1)+y(n+1)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));
%再在指令窗口执行
ex5_12(inline('1/sqrt(2*pi)*exp(-x.^2/2)'),0,1,50)
%Exercise 13
fun=inline('5400*v./(8.276*v.^2+2000)','v');
quadl(fun,15,30)
%Exercise 14
重心不超过凳边沿。1/2, 2/3, 3/4, ...,n/(n+1)
%Exercise15
利润函数fun=inline('(p-c0+k*log(M*exp(-a*p)))*M*exp(-a*p)','p');
求p使fun最大
%Exercise 16
clear; x=-3/4:0.01:3/4;
y=(3/4+x)*2.*sqrt(1-16/9.*x.^2)*9.8;
P=trapz(x,y) %单位:千牛
%Exercise 17
clear; close;
fplot('17-t^(2/3)-5-2*t^(2/3)',[0,20]); grid;
t=fzero('17-x^(2/3)-5-2*x^(2/3)',7)
t=0:0.1:8; y=17-t.^(2/3)-5-2*t.^(2/3);
trapz(t,y)-20 %单位:百万元
%Exercise 18
%曲面面积计算
%Excercise 1(1)
fun=inline('x+y','x','y');
[t,y]=ode45(fun,[0 1 2 3],1) %注意由于初值为y(0)=1,[0 1 2 3]中0不可缺
%Excercise 1(3)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=0.01*y(2)^2-2*y(1)+sin(t)
%运行下列指令
clear;close;
fun=@(t,y)[y(2);0.01*y(2)^2-2*y(1)+sin(t)];
[t,y]=ode45(fun,[0 5],[0;1]);
plot(t,y(:,1))
%Excercise 1(5)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=-mu*(y(1)^2-1)*y(2)-y(1)
%运行下列指令,注意参数mu的处理
clear;close;
fun=@(t,y,mu)[y(2);-mu*(y(1)^2-1)*y(2)-y(1)];
[t,y]=ode45(fun,[0 20],[2;0],[],1);
plot(y(:,1),y(:,2));hold on;
[t,y]=ode45(fun,[0 20],[2;0],[],2);
plot(y(:,1),y(:,2),'r');hold off;
%Excercise 2
roots([1 10 54 132 137 50])
%通解A1*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-t)+A5*t*exp(-t) %Excercise 3
dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');
[x,y]=ode45(dfun,[0,50],[1;-1]);length(x)
%所用节点很多
[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x)
%所用节点很少
%Excercise 4
clear;
dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3;
x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3]','t','x');
[t,x]=ode45(dfun,[0 24],[1.2; 0; 0; -1.04935371]);
plot(x(:,1),x(:,3));
%Excercise 5
%方程y'=2x+y^2,y(0)=0
clear;close;
fun=inline('2*x+y^2','x','y');
[x,y]=ode45(fun,[0 1.57],0); %x的上界再增加,解会
plot(x,y)
%Excercise 6
clear;close;
fun=@(t,x,a,b)a*x+b;
[t,x]=ode45(fun,[0 10],0.1,[],1,1);
subplot(2,4,1);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,1);
subplot(2,4,2);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],1,-1);
subplot(2,4,3);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,-1);
subplot(2,4,4);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,1);
subplot(2,4,5);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,1);
subplot(2,4,6);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,-1);
subplot(2,4,7);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,-1);
subplot(2,4,8);plot(t,x)
%Excercise 7
%微分方程 T'=k(c-T),T(0)=20
dsolve('DT=k*(c-T)','T(0)=20','t')
%得c+exp(-k*t)*(-c+20)
%利用T(10)=25.2, T(20)=28.32拟合(或者解非线性方程)
fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t')
lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32])
%解得户外温度c=33,比例系数k=0.05.
%Excercise 8
%微分方程 x'/x=0.5*(1-x),x(0)=0.05
fun=inline('0.5*(1-x)*x','t','x');
[t,x]=ode45(fun,[0 10],0.05);
plot(t,x)
id=min(find(x>0.8));
t(id)
%Excercise 9
%微分方程组 V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t)
%答案(1)exp(20);(2)0.353;(3)30;(4)451,0.4,9.6
%Chapter 7
%Exercise 1
syms ph th;
a=sin(ph)*cos(th)-cos(ph)*sin(th)-sin(ph-th);
simple(a)
%化简后差的结果为0
%Exercise 2
syms x;s=x^4-5*x^3+5*x^2+5*x-6;
factor(s)
%Exercise 3
syms a;A=[1 2;2 a];
iA=inv(A),[v,d]=eig(A)
%Exercise 4
syms x y;
limit((3^x+9^x)^(1/x),x,inf)
s1=limit(x*y/(sqrt(x*y+1)-1),x,0);s2=limit(s1,y,0)
%Exercise 5
syms k n x;s1=symsum(k^2,k,1,n);s1=simple(s1)
s2=symsum(k^(-2),k,1,inf);s2=simple(s2)
s3=symsum(1/(2*n+1)/(2*x+1)^(2*n+1),n,0,inf);s3=simple(s3) %Exercise 6
syms x y z;s=sin(x^2*y*z);
s=diff(s,x,2);
s=diff(s,y,1);
s=subs(s,{x,y,z},{1,1,3})
%Exercise 7
syms x;s=log(x+sqrt(1+x^2));taylor(s,8,0,x)
%Exercise 8 (以第四章习题9为例)
先用符号运算求偏导数
syms x y;f=y^3/9+3*x^2*y+9*x^2+y^2+x*y+9;
fx=diff(f,x),fy=diff(f,y)
根据计算结果得方程组. 求解方程组
[sx,sy]=solve(6*x*y+18*x+y,1/3*y^2+3*x^2+2*y+x,x,y)
得四个解(0,0),(-1/3,-6),(-7/6,-7/2),(5/6,-5/2).计算Hesse矩阵 fh2=[diff(fx,x),diff(fx,y);diff(fy,x),diff(fy,y)]
计算
eig(subs(fh2,[x,y],[0,0]))
得知正定,所以是极小值点.极小值用
subs(f,[x,y],[0,0])
求得。同理可得(-1/3,-6)为极大值点,其它两个为鞍点。 %Exercise 9(以第一小题为例)
syms y;f=exp(2*y)/(exp(y)+2);
fi=int(f,y)
s=simple(diff(fi)-f)
%Exercise 10
syms x y;f=(x-y)^3*sin(x+2*y);Ix=simple(int(f,y,-x,x))
%Exercise 11(3)
syms x;f=x*log(x^4)*asin(1/x);Ix=int(f,x,1,3);vpa(Ix)
%Exercise 12
%1(3)
syms x;solve(5*x^23-6*x^7+8*x^6-5*x^2)
%6
syms a b;s=solve(a-0.7*sin(a)-0.2*cos(b),b-0.7*cos(a)+0.2*sin(b)); s.a,s.b
%Exercise 13
%1(3)
dsolve('D2y-0.01*Dy^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')(解不出) %1(4)
dsolve('2*D2x-5*Dx+3*x=45*exp(2*t)','x(0)=2','Dx(0)=1','t') %Exercise 14
%6(ii)
ezplot('x^2/4+y^2/9=1')
%6(vi)
ezmesh('2*sin(ph)*cos(th)','2*sin(ph)*sin(th)','2*cos(ph)',[0 pi/2 0 2*pi])
范文五:matlab数学实验答案(胡良剑)
%Page20,ex1
(5) 等于[exp(1),exp(2);exp(3),exp(4)] (7) 3=1*3, 8=2*4
(8) a为各列最小值,b为最小值所在的行号
(10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture
(11) 答案表明:编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)
(12) 答案表明:编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不
等式(40>=10)
%Page20, ex2
(1)a, b, c的值尽管都是1, 但数据类型分别为数值,字符, 逻辑, 注意a与c相等, 但
他们不等于b
(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码 %Page20,ex3
>> r=2;p=0.5;n=12;
>> T=log(r)/n/log(1+0.01*p) T =
11.5813
%Page20,ex4
>> x=-2:0.05:2;f=x.^4-2.^x; >> [fmin,min_index]=min(f) fmin =
-1.3907 %最小值
min_index =
54 %最小值点编址
>> x(min_index)
ans =
0.6500 %最小值点
>> [f1,x1_index]=min(abs(f)) %求近似根--绝对值最小的点 f1 =
0.0328
x1_index =
24
>> x(x1_index)
ans =
-0.8500
>> x(x1_index)=[];f=x.^4-2.^x; %删去绝对值最小的点以求函数绝对值次小的点 >> [f2,x2_index]=min(abs(f)) %求另一近似根--函数绝对值次小的点 f2 =
0.0630
x2_index =
65
>> x(x2_index)
ans =
1.2500
%Page20,ex5
>> z=magic(10)
z =
92 99 1 8 15 67 74 51 58 40 98 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 56 63 70 47 85 87 19 21 3 60 62 69 71 28 86 93 25 2 9 61 68 75 52 34 17 24 76 83 90 42 49 26 33 65 23 5 82 89 91 48 30 32 39 66 79 6 13 95 97 29 31 38 45 72 10 12 94 96 78 35 37 44 46 53 11 18 100 77 84 36 43 50 27 59 >> sum(z)
ans =
505 505 505 505 505 505 505 505 505 505
>> sum(diag(z))
ans =
505
>> z(:,2)/sqrt(3)
ans =
57.1577
46.1880
46.7654
50.2295
53.6936
13.8564
2.8868
3.4641
6.9282
10.3923
>> z(8,:)=z(8,:)+z(3,:) z =
92 99 1 8 15 67 74 51 58 40 98 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 56 63 70 47 85 87 19 21 3 60 62 69 71 28 86 93 25 2 9 61 68 75 52 34 17 24 76 83 90 42 49 26 33 65 23 5 82 89 91 48 30 32 39 66
83 87 101 115 119 83 87 101 115 119
10 12 94 96 78 35 37 44 46 53
11 18 100 77 84 36 43 50 27 59
%Page 40 ex1
先在编辑器窗口写下列M函数,保存为eg2_1.m function [xbar,s]=ex2_1(x) n=length(x);
xbar=sum(x)/n;
s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));
例如
>>x=[81 70 65 51 76 66 90 87 61 77];
>>[xbar,s]=ex2_1(x) xbar =
72.4000
s =
12.1124
%Page 40 ex2
s=log(1);n=0;
while s<=100>=100>
n=n+1;
s=s+log(1+n);
end
m=n
计算结果m=37
%Page 40 ex3
clear;
F(1)=1;F(2)=1;k=2;x=0; e=1e-8; a=(1+sqrt(5))/2; while abs(x-a)>e
k=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);
end
a,x,k
计算至k=21可满足精度
%Page 40 ex4
clear;tic;s=0;
for i=1:1000000
s=s+sqrt(3)/2^i;
end
s,toc
tic;s=0;i=1;
while i<=1000000>=1000000>
s=s+sqrt(3)/2^i;i=i+1;
end
s,toc
tic;s=0;
i=1:1000000;
s=sqrt(3)*sum(1./2.^i);
s,toc
%Page 40 ex5
t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16]; plot(t,c)
%Page 40 ex6
%(1)
x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)
y=inline('x^2*sin(x^2-x-2)');fplot(y,[-2 2])
%(2)参数方法
t=linspace(0,2*pi,100);
x=2*cos(t);y=3*sin(t); plot(x,y) %(3)
x=-3:0.1:3;y=x;
[x,y]=meshgrid(x,y);
z=x.^2+y.^2;
surf(x,y,z)
%(4)
x=-3:0.1:3;y=-3:0.1:13;
[x,y]=meshgrid(x,y);
z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6; surf(x,y,z)
%(5)
t=0:0.01:2*pi;
x=sin(t);y=cos(t);z=cos(2*t); plot3(x,y,z)
%(6)
theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20);
[theta,fai]=meshgrid(theta,fai); x=2*sin(fai).*cos(theta);
y=2*sin(fai).*sin(theta);z=2*cos(fai); surf(x,y,z)
%(7)
x=linspace(0,pi,100);
y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x); plot(x,y1,x,y2,x,y3)
%page41, ex7
x=-1.5:0.05:1.5;
y=1.1*(x>1.1)+x.*(x<=1.1).*(x>=-1.1)-1.1*(x<-1.1); plot(x,y)="">-1.1);>
%page41,ex8
分别使用which trapz, type trapz, dir C:\MATLAB7\toolbox\matlab\datafun\
%page41,ex9
clear;close;
x=-2:0.1:2;y=x;
[x,y]=meshgrid(x,y);
a=0.5457;b=0.7575;
p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1); p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y<=1); p="">=1);><=-1); mesh(x,y,p)="">=-1);>
%page41, ex10
lookfor lyapunov
help lyap
>> A=[1 2 3;4 5 6;7 8 0];C=[2 -5 -22;-5 -24 -56;-22 -56 -16];
>> X=lyap(A,C)
X =
1.0000 -1.0000 -0.0000
-1.0000 2.0000 1.0000
-0.0000 1.0000 7.0000
%Chapter 3
%Exercise 1
>> a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\b
ans =
0.5000 0.5000 1.0000 ans =
2 2 1
ans =
0.6552 %一元方程组x[2,4,3]=[1,2,3]的近似解 ans =
0 0 0
0 0 0
0.6667 1.3333 1.0000 %矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解
%Exercise 2(1)
>> A=[4 1 -1;3 2 -6;1 -5 3];b=[9;-2;1];
>> rank(A), rank([A,b]) %[A,b]为增广矩阵 ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
2.3830
1.4894
2.0213
%Exercise 2(2)
>> A=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];
>> rank(A), rank([A,b]) ans =
3
ans =
3 %可见方程组唯一解
>> x=A\b
x =
-0.4706
-0.2941
0
%Exercise 2(3)
>> A=[4 1;3 2;1 -5];b=[1;1;1];
>> rank(A), rank([A,b])
ans =
2
ans =
3 %可见方程组无解
>> x=A\b
x =
0.3311
-0.1219 %最小二乘近似解
%Exercise 2(4)
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1 2 3]';%注意b的写法
>> rank(a),rank([a,b]) ans =
3
ans =
3 %rank(a)==rank([a,b])<4说明有无穷多解>> a\b
ans =
1
0
1
0 %一个特解
%Exercise 3
>> a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';
>> x=null(a),x0=a\b x =
-0.6255
0.6255
-0.2085
0.4170
x0 =
1
0
1
0
%通解kx+x0
%Exercise 4
>> x0=[0.2 0.8]';a=[0.99 0.05;0.01 0.95];
>> x1=a*x, x2=a^2*x, x10=a^10*x
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> x0=[0.8 0.2]';
>> x=x0;for i=1:1000,x=a*x;end,x
x =
0.8333
0.1667
>> [v,e]=eig(a)
v =
0.9806 -0.7071
0.1961 0.7071
e =
1.0000 0
0 0.9400
>> v(:,1)./x
ans =
1.1767
1.1767 %成比例,说明x是最大特征值对应的特征向量 %Exercise 5
%用到公式(3.11)(3.12)
>> B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25 5 20]';
>> C=B/diag(x)
C =
0.2400 0.4000 0.0500 0.0900 0.2000 0.0100 0.1200 0.0400 0.0900 >> A=eye(3,3)-C
A =
0.7600 -0.4000 -0.0500 -0.0900 0.8000 -0.0100 -0.1200 -0.0400 0.9100 >> D=[17 17 17]';x=A\D x =
37.5696
25.7862
24.7690
%Exercise 6(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a)
ans =
-94
ans =
0.2553 -0.0213 0.0426 0.1596 -0.1383 -0.2234 0.1809 -0.2234 -0.0532 v =
0.0185 -0.9009 -0.3066 -0.7693 -0.1240 -0.7248 -0.6386 -0.4158 0.6170 d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
%Exercise 6(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];det(a),inv(a),[v,d]=eig(a)
ans =
1
ans =
2.0000 -2.0000 1.0000 1.0000 -1.0000 1.0000 2.0000 -3.0000 2.0000 v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774 -0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0 0 0 1.0000 - 0.0000i %Exercise 6(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> det(A),inv(A), [v,d]=eig(A)
ans =
1
ans =
68.0000 -41.0000 -17.0000 10.0000 -41.0000 25.0000 10.0000 -6.0000 -17.0000 10.0000 5.0000 -3.0000 10.0000 -6.0000 -3.0000 2.0000 v =
0.8304 0.0933 0.3963 0.3803 -0.5016 -0.3017 0.6149 0.5286 -0.2086 0.7603 -0.2716 0.5520 0.1237 -0.5676 -0.6254 0.5209 d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
%Exercise 6(4)(以n=5为例)
%关键是矩阵的定义
%方法一(三个for)
n=5;
for i=1:n, a(i,i)=5;end
for i=1:(n-1),a(i,i+1)=6;end for i=1:(n-1),a(i+1,i)=1;end a
%方法二(一个for)
n=5;a=zeros(n,n);
a(1,1:2)=[5 6];
for i=2:(n-1),a(i,[i-1,i,i+1])=[1 5 6];end
a(n,[n-1 n])=[1 5];
a
%方法三(不用for)
n=5;a=diag(5*ones(n,1)); b=diag(6*ones(n-1,1));
c=diag(ones(n-1,1));
a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)]
%下列计算
>> det(a)
ans =
665
>> inv(a)
ans =
0.3173 -0.5865 1.0286 -1.6241 1.9489 -0.0977 0.4887 -0.8571 1.3534 -1.6241
0.0286 -0.1429 0.5429 -0.8571 1.0286 -0.0075 0.0376 -0.1429 0.4887 -0.5865 0.0015 -0.0075 0.0286 -0.0977 0.3173 >> [v,d]=eig(a)
v =
-0.7843 -0.7843 -0.9237 0.9860 -0.9237 0.5546 -0.5546 -0.3771 -0.0000 0.3771 -0.2614 -0.2614 0.0000 -0.1643 0.0000 0.0924 -0.0924 0.0628 -0.0000 -0.0628 -0.0218 -0.0218 0.0257 0.0274 0.0257 d =
0.7574 0 0 0 0
0 9.2426 0 0 0
0 0 7.4495 0 0
0 0 0 5.0000 0
0 0 0 0 2.5505
%Exercise 7(1)
>> a=[4 1 -1;3 2 -6;1 -5 3];[v,d]=eig(a) v =
0.0185 -0.9009 -0.3066
-0.7693 -0.1240 -0.7248
-0.6386 -0.4158 0.6170
d =
-3.0527 0 0
0 3.6760 0
0 0 8.3766
>> det(v)
ans =
-0.9255 %v行列式正常, 特征向量线性相关,可对角化
>> inv(v)*a*v %验算
ans =
-3.0527 0.0000 -0.0000
0.0000 3.6760 -0.0000
-0.0000 -0.0000 8.3766
>> [v2,d2]=jordan(a) %也可用jordan
v2 =
0.0798 0.0076 0.9127
0.1886 -0.3141 0.1256
-0.1605 -0.2607 0.4213 %特征向量不同
d2 =
8.3766 0 0
0 -3.0527 - 0.0000i 0
0 0 3.6760 + 0.0000i >> v2\a*v2
ans =
8.3766 0 0.0000
0.0000 -3.0527 0.0000 0.0000 0.0000 3.6760 >> v(:,1)./v2(:,2) %对应相同特征值的特征向量成比例 ans =
2.4491
2.4491
2.4491
%Exercise 7(2)
>> a=[1 1 -1;0 2 -1;-1 2 0];[v,d]=eig(a)
v =
-0.5773 0.5774 + 0.0000i 0.5774 - 0.0000i
-0.5773 0.5774 0.5774 -0.5774 0.5773 - 0.0000i 0.5773 + 0.0000i
d =
1.0000 0 0
0 1.0000 + 0.0000i 0 0 0 1.0000 - 0.0000i >> det(v)
ans =
-5.0566e-028 -5.1918e-017i %v的行列式接近0, 特征向量线性相关,不可对角化
>> [v,d]=jordan(a) v =
1 0 1
1 0 0
1 -1 0
d =
1 1 0
0 1 1
0 0 1 %jordan标准形不是对角的,所以不可对角化 %Exercise 7(3)
>> A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
A =
5 7 6 5
7 10 8 7
6 8 10 9
5 7 9 10
>> [v,d]=eig(A)
v =
0.8304 0.0933 0.3963 0.3803 -0.5016 -0.3017 0.6149 0.5286 -0.2086 0.7603 -0.2716 0.5520 0.1237 -0.5676 -0.6254 0.5209 d =
0.0102 0 0 0
0 0.8431 0 0
0 0 3.8581 0
0 0 0 30.2887
>> inv(v)*A*v
ans =
0.0102 0.0000 -0.0000 0.0000 0.0000 0.8431 -0.0000 -0.0000 -0.0000 0.0000 3.8581 -0.0000 -0.0000 -0.0000 0 30.2887 %本题用jordan不行, 原因未知
%Exercise 7(4)参考6(4)和7(1), 略
%Exercise 8 只有(3)对称, 且特征值全部大于零, 所以是正定矩阵. %Exercise 9(1)
>> a=[4 -3 1 3;2 -1 3 5;1 -1 -1 -1;3 -2 3 4;7 -6 -7 0]
>> rank(a)
ans =
3
>> rank(a(1:3,:))
ans =
2
>> rank(a([1 2 4],:)) %1,2,4行为最大无关组
ans =
3
>> b=a([1 2 4],:)';c=a([3 5],:)'; >> b\c %线性表示的系数
ans =
0.5000 5.0000
-0.5000 1.0000
0 -5.0000
%Exercise 10
>> a=[1 -2 2;-2 -2 4;2 4 -2] >> [v,d]=eig(a)
v =
0.3333 0.9339 -0.1293 0.6667 -0.3304 -0.6681 -0.6667 0.1365 -0.7327 d =
-7.0000 0 0
0 2.0000 0
0 0 2.0000
>> v'*v
ans =
1.0000 0.0000 0.0000 0.0000 1.0000 0
0.0000 0 1.0000 %v确实是正交矩阵
%Exercise 11
%设经过6个电阻的电流分别为i1, ..., i6. 列方程组如下 %20-2i1=a; 5-3i2=c; a-3i3=c; a-4i4=b; c-5i5=b; b-3i6=0;
%i1=i3+i4;i5=i2+i3;i6=i4+i5; %计算如下
>> A=[1 0 0 2 0 0 0 0 0;0 0 1 0 3 0 0 0 0;1 0 -1 0 0 -3 0 0 0; 1 -1 0 0 0 0 -4 0 0;
0 -1 1 0 0 0 0 -5 0;0 1 0 0 0 0 0 0 -3; 0 0 0 1 0 -1 -1 0 0;0 0 0 0 -1 -1 0 1 0;
0 0 0 0 0 0 -1 -1 1]; >>b=[20 5 0 0 0 0 0 0 0]'; A\b ans =
13.3453
6.4401
8.5420
3.3274
-1.1807
1.6011
1.7263
0.4204
2.1467
%Exercise 12
>> A=[1 2 3;4 5 6;7 8 0]; >> left=sum(eig(A)), right=sum(trace(A))
left =
6.0000
right =
6
>> left=prod(eig(A)), right=det(A) %原题有错, (-1)^n应删去 left =
27.0000
right =
27
>> fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3))
fA =
1.0e-012 *
0.0853 0.1421 0.0284 0.1421 0.1421 0
-0.0568 -0.1137 0.1705 >> norm(fA) %f(A)范数接近0
ans =
2.9536e-013
%Exercise 1(1)
roots([1 1 1])
%Exercise 1(2)
roots([3 0 -4 0 2 -1])
%Exercise 1(3)
p=zeros(1,24);
p([1 17 18 22])=[5 -6 8 -5]; roots(p)
%Exercise 1(4)
p1=[2 3];
p2=conv(p1, p1);
p3=conv(p1, p2);
p3(end)=p3(end)-4; %原p3最后一个分量-4
roots(p3)
%Exercise 2
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
fzero(fun,2)
%Exercise 3
fun=inline('x^4-2^x');
fplot(fun,[-2 2]);grid on; fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)
%Exercise 4
fun=inline('x*sin(1/x)','x'); fplot(fun, [-0.1 0.1]);
x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0.5)*0.01);end;
x=[x,-x]
%Exercise 5
fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^
2-16*x(3)^2]','x');
[a,b,c]=fsolve(fun,[0 0 0]) %Exercise 6
fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];
[a,b,c]=fsolve(fun,[0.5 0.5]) %Exercise 7
clear; close; t=0:pi/100:2*pi; x1=2+sqrt(5)*cos(t); y1=3-2*x1+sqrt(5)*sin(t);
x2=3+sqrt(2)*cos(t); y2=6*sin(t); plot(x1,y1,x2,y2); grid on; %作图发现4个解的大致位置,然后分别求解
y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])
y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])
y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])
y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])
%Exercise 8(1)
clear;
fun=inline('x.^2.*sin(x.^2-x-2)'); fplot(fun,[-2 2]);grid on; %作图观察
x(1)=-2;
x(3)=fminbnd(fun,-1,-0.5);
x(5)=fminbnd(fun,1,2);
fun2=inline('-x.^2.*sin(x.^2-x-2)'); x(2)=fminbnd(fun2,-2,-1);
x(4)=fminbnd(fun2,-0.5,0.5);
x(6)=2
feval(fun,x)
%答案: 以上x(1)(3)(5)是局部极小,x(2)(4)(6)是局部极大,从最后一句知道x(1)全局最小, x(2)最大。
%Exercise 8(2)
clear;
fun=inline('3*x.^5-20*x.^3+10'); fplot(fun,[-3 3]);grid on;%作图观察
x(1)=-3;
x(3)=fminsearch(fun,2.5);
fun2=inline('-(3*x.^5-20*x.^3+10)'); x(2)=fminsearch(fun2,-2.5);
x(4)=3;
feval(fun,x)
%Exercise 8(3)
fun=inline('abs(x^3-x^2-x-2)'); fplot(fun,[0 3]);grid on;%作图观察
fminbnd(fun,1.5,2.5)
fun2=inline('-abs(x^3-x^2-x-2)'); fminbnd(fun2,0.5,1.5)
%Exercise 9
close;
x=-2:0.1:1;y=-7:0.1:1;
[x,y]=meshgrid(x,y);
z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9; mesh(x,y,z);grid on;%作图观察
fun=inline('x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9');
x=fminsearch(fun,[0 0])%求极小值
fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)');
x=fminsearch(fun2,[0 -5])%求极大值
%Exercise 10
clear;t=0:24;
c=[15 14 14 14 14 15 16 18 20 22 23 25 28 ...
31 32 31 29 27 25 24 22 20 18 17 16]; p2=polyfit(t,c,2)
p3=polyfit(t,c,3)
fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t');
a=lsqcurvefit(fun,[0 0],t,c)%初值可以试探
f=feval(fun, a,t)
norm(f-c)%拟合效果
plot(t,c,t,f) %作图检验
fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改f(x)+20 b=lsqcurvefit(fun2,[0 0],t,c) figure
f2=feval(fun2, b,t)
norm(f2-c)%拟合效果
plot(t,c,t,f2) %作图检验
%Exercise 11
fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)');
x=fzero(fun, 0, 1)
%Exercise 12
r=5.04/12/100;N=20*12;
x=7500*180 %房屋总价格
y=x*0.3 %首付款额
x0=x-y%贷款总额
a=(1+r)^N*r*x0/((1+r)^N-1)%月付还款额
r1=4.05/12/100;x1=10*10000;%公积金贷款
a1=(1+r1)^N*r1*x1/((1+r1)^N-1) x2=x0-x1%商业贷款
a2=(1+r)^N*r*x2/((1+r)^N-1)
a=a1+a2
%Exercise 13
%列方程th*R^2+(pi-2*th)*r^2-R*r*sin(th)=pi*r^2/2
%化简得sin(2*th)-2*th*cos(2*th)=pi/2 %以下Matlab计算
clear;fun= inline('sin(2*th)-2*th*cos(2*th)-pi/2','th')
th=fsolve(fun,pi/4)
R=20*cos(th)
%Exercise 14
%先在Editor窗口写M函数保存
function x=secant(fname,x0,x1,e) while abs(x0-x1)>e,
x=x1-(x1-x0)*feval(fname,x1)/(feval(fname,x1)-feval(fname,x0));
x0=x1;x1=x;
end
%再在指令窗口
fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');
secant(fun,1,2,1e-8)
%Exercise 15
%作系数为a,初值为xo,从第m步到第n步迭代过程的M函数:
function f=ex4_15fun(a,x0,m,n) x(1)=x0; y(1)=a*x(1)+1;x(2)=y(1); if m<2, plot([x(1),x(1),x(2)],[0,y(1),y(1)]);hold="" on;="" end="">2,>
for i=2:n
y(i)=a*x(i)+1; x(i+1)=y(i);
if i>m, plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]); end
end
hold off;
%M脚本文件
subplot(2,2,1);ex4_15fun(0.9,1,1,20); subplot(2,2,2);ex4_15fun(-0.9,1,1,20); subplot(2,2,3);ex4_15fun(1.1,1,1,20); subplot(2,2,4);ex4_15fun(-1.1,1,1,20); %Exercise 16
%设夹角t, 问题转化为 min f=5/sin(t)+10/cos(t) %取初始值pi/4, 计算如下
fun=@(t)5/sin(t)+10/cos(t);
[t,f]=fminsearch(fun, pi/4)
t =
0.6709
f =
20.8097
%Exercise 17
%提示:x(k+2)=f(x(k))=a^2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k)))
%计算平衡点x
%|f'(x)|<1则稳定>1则稳定>
%Exercise 18
%先写M文件
function f=ex4_18(a,x0,n)
x=zeros(1,n);y=x;
x(1)=x0;
y(1)=a*x(1)+1;
x(2)=y(1);
plot([x(1),x(1),x(2)],[0,y(1),y(1)],'r'); hold on;
for i=2:n
y(i)=a*x(i)+1;
x(i+1)=y(i);
plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]) end
hold off;
%再执行指令
>> ex4_18(0.9,1,20)
>> ex4_18(-0.9,1,20)
>> ex4_18(1.1,1,20)
>> ex4_18(-1.1,1,20)
%Exercise 19
clear; close; x(1)=0; y(1)=0;
for k=1:3000
x(k+1)=1+y(k)-1.4*x(k)^2; y(k+1)=0.3*x(k); end
plot(x(1000:1500),y(1000:1500),'+g');hold on plot(x(1501:2000),y(1501:2000),'.b'); plot(x(2001:2500),y(2001:2500),'*y'); plot(x(2501:3001),y(2501:3001),'.r');
%Exercise 1
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0]; trapz(x,y)
%Exercise 2
x=[0 4 10 12 15 22 28 34 40];
y=[0 1 3 6 8 9 5 3 0]; diff(y)./diff(x)
%Exercise 3
xa=-1:0.1:1;ya=0:0.1:2; [x,y]=meshgrid(xa,ya); z=x.*exp(-x.^2 -y.^3); [px,py] = gradient(z,xa,ya); px
%Exercise 4
t=0:0.01:1.5;
x=log(cos(t));
y=cos(t)-t.*sin(t); dydx=gradient(y,x) [x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点 dydx(id)
%Exercise 5(2)
fun=inline('exp(2*x).*cos(x).^3');
quadl(fun,0,2*pi)
或用trapz
x=linspace(0,2*pi,100); y=exp(2*x).*cos(x).^3; trapz(x,y)
%Exercise 5(3)
fun=@(x)x.*log(x.^4).*asin(1./x.^2);
quadl(fun,1,3)
或用trapz
x=1:0.01:3;
y=feval(fun,x);
trapz(x,y)
%Exercise 5(4)
fun=@(x)sin(x)./x;
quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义,用很小的1e-10代替
%Exercise 5(5)
%参考Exercise 5(4)
%Exercise 5(6)
fun=inline('sqrt(1+r.^2.*sin(th))','r','th');
dblquad(fun,0,1,0,2*pi) %Exercise 5(7)
首先建立84页函数dblquad2
clear;
fun=@(x,y)1+x+y.^2;
clo=@(x)-sqrt(2*x-x.^2); dup=@(x)sqrt(2*x-x.^2); dblquad2(fun,0,2,clo,dhi,100) %Exercise 6
t=linspace(0,2*pi,100); x=2*cos(t);y=3*sin(t); dx=gradient(x,t);dy=gradient(y,t);
f=sqrt(dx.^2+dy.^2);
trapz(t,f)
%Exercise 7
xa=-1:0.1:1;ya=0:0.1:2; [x,y]=meshgrid(xa,ya); z=x.*exp(x.^2+y.^2);
[zx,zy]=gradient(z,xa,ya); f=sqrt(1+zx.^2+zy.^2); s=0;
for i=2:length(xa)
for j=2:length(ya)
s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4;
end
end
s
%Exercise 8
funl=inline('-(-x).^0.2.*cos(x)');
funr=inline('x.^0.2.*cos(x)');
quadl(funl,-1,0)+quadl(funr,0,1) %Exercise 9 (以I32为例)
fun=@(x)abs(sin(x));
h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y) h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致,结果失真
q1=quad(fun,0,32*pi)
q2=quadl(fun,0,32*pi)
%Exercise 10(2)
%先在程序编辑器,写下列函数,保存为ex5_10_2f
function d=ex5_10_2f(fname,a,h0,e) h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);
d0=d+2*e;
while abs(d-d0)>e
d0=d;h0=h;h=h0/2;
d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);
end
%再在指令窗口执行
fun=inline('x.^2*sin(x.^2-x-2)','x'); d=ex5_10_2f(fun,1.4,0.1,1e-3)
%Exercise 11
%提示:f上升时,f'>0;f下降时,f'<0; f极值,="" f'="0.">0;>
%Exercise 12
在程序编辑器,写下列函数,保存为ex5_12f
function I=ex5_12(fname,a,b,n)
x=linspace(a,b,n+1);
y=feval(fname,x);
I=(b-a)/n/3*(y(1)+y(n+1)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));
%再在指令窗口执行
ex5_12(inline('1/sqrt(2*pi)*exp(-x.^2/2)'),0,1,50)
%Exercise 13
fun=inline('5400*v./(8.276*v.^2+2000)','v'); quadl(fun,15,30)
%Exercise 14
重心不超过凳边沿。1/2, 2/3, 3/4, ...,n/(n+1)
%Exercise15
利润函数fun=inline('(p-c0+k*log(M*exp(-a*p)))*M*exp(-a*p)','p');
求p使fun最大
%Exercise 16
clear; x=-3/4:0.01:3/4;
y=(3/4+x)*2.*sqrt(1-16/9.*x.^2)*9.8; P=trapz(x,y) %单位:千牛
%Exercise 17
clear; close;
fplot('17-t^(2/3)-5-2*t^(2/3)',[0,20]); grid; t=fzero('17-x^(2/3)-5-2*x^(2/3)',7) t=0:0.1:8; y=17-t.^(2/3)-5-2*t.^(2/3); trapz(t,y)-20 %单位:百万元
%Exercise 18
%曲面面积计算
%Excercise 1(1)
fun=inline('x+y','x','y');
[t,y]=ode45(fun,[0 1 2 3],1) %注意由于初值为y(0)=1,[0 1 2 3]中0不可缺
%Excercise 1(3)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=0.01*y(2)^2-2*y(1)+sin(t) %运行下列指令
clear;close;
fun=@(t,y)[y(2);0.01*y(2)^2-2*y(1)+sin(t)]; [t,y]=ode45(fun,[0 5],[0;1]);
plot(t,y(:,1))
%Excercise 1(5)
%令y(1)=y,y(2)=y',化为方程组
%y(1)'=y(2),y(2)'=-mu*(y(1)^2-1)*y(2)-y(1) %运行下列指令,注意参数mu的处理
clear;close;
fun=@(t,y,mu)[y(2);-mu*(y(1)^2-1)*y(2)-y(1)]; [t,y]=ode45(fun,[0 20],[2;0],[],1); plot(y(:,1),y(:,2));hold on;
[t,y]=ode45(fun,[0 20],[2;0],[],2); plot(y(:,1),y(:,2),'r');hold off;
%Excercise 2
roots([1 10 54 132 137 50])
%通解A1*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-t)+A5*t*exp(-t)
%Excercise 3
dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');
[x,y]=ode45(dfun,[0,50],[1;-1]);length(x) %所用节点很多
[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x) %所用节点很少
%Excercise 4
clear;
dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(
1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3;
x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1
)+1-1/82.45)^2+x(3)^2))^3]','t','x'); [t,x]=ode45(dfun,[0 24],[1.2; 0; 0; -1.04935371]); plot(x(:,1),x(:,3));
%Excercise 5
%方程y'=2x+y^2,y(0)=0
clear;close;
fun=inline('2*x+y^2','x','y');
[x,y]=ode45(fun,[0 1.57],0); %x的上界再增加,解会"爆炸"
plot(x,y)
%Excercise 6
clear;close;
fun=@(t,x,a,b)a*x+b;
[t,x]=ode45(fun,[0 10],0.1,[],1,1); subplot(2,4,1);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,1); subplot(2,4,2);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],1,-1); subplot(2,4,3);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],1,-1); subplot(2,4,4);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,1); subplot(2,4,5);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,1); subplot(2,4,6);plot(t,x)
[t,x]=ode45(fun,[0 10],0.1,[],-1,-1); subplot(2,4,7);plot(t,x)
[t,x]=ode45(fun,[0 10],-0.1,[],-1,-1); subplot(2,4,8);plot(t,x)
%Excercise 7
%微分方程 T'=k(c-T),T(0)=20
dsolve('DT=k*(c-T)','T(0)=20','t') %得c+exp(-k*t)*(-c+20)
%利用T(10)=25.2, T(20)=28.32拟合(或者解非线性方程)
fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t') lsqcurvefit(fun,[30 1],[10 20],[25.2 28.32]) %解得户外温度c=33,比例系数k=0.05.
%Excercise 8
%微分方程 x'/x=0.5*(1-x),x(0)=0.05
fun=inline('0.5*(1-x)*x','t','x'); [t,x]=ode45(fun,[0 10],0.05);
plot(t,x)
id=min(find(x>0.8));
t(id)
%Excercise 9
%微分方程组 V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t) %答案(1)exp(20);(2)0.353;(3)30;(4)451,0.4,9.6
%Chapter 7
%Exercise 1
syms ph th;
a=sin(ph)*cos(th)-cos(ph)*sin(th)-sin(ph-th); simple(a)
%化简后差的结果为0
%Exercise 2
syms x;s=x^4-5*x^3+5*x^2+5*x-6;
factor(s)
%Exercise 3
syms a;A=[1 2;2 a];
iA=inv(A),[v,d]=eig(A)
%Exercise 4
syms x y;
limit((3^x+9^x)^(1/x),x,inf)
s1=limit(x*y/(sqrt(x*y+1)-1),x,0);s2=limit(s1,y,0) %Exercise 5
syms k n x;s1=symsum(k^2,k,1,n);s1=simple(s1) s2=symsum(k^(-2),k,1,inf);s2=simple(s2) s3=symsum(1/(2*n+1)/(2*x+1)^(2*n+1),n,0,inf);s3=simple(s3)
%Exercise 6
syms x y z;s=sin(x^2*y*z);
s=diff(s,x,2);
s=diff(s,y,1);
s=subs(s,{x,y,z},{1,1,3})
%Exercise 7
syms x;s=log(x+sqrt(1+x^2));taylor(s,8,0,x) %Exercise 8 (以第四章习题9为例)
先用符号运算求偏导数
syms x y;f=y^3/9+3*x^2*y+9*x^2+y^2+x*y+9; fx=diff(f,x),fy=diff(f,y)
根据计算结果得方程组. 求解方程组
[sx,sy]=solve(6*x*y+18*x+y,1/3*y^2+3*x^2+2*y+x,x,y)
得四个解(0,0),(-1/3,-6),(-7/6,-7/2),(5/6,-5/2).计算Hesse矩阵
fh2=[diff(fx,x),diff(fx,y);diff(fy,x),diff(fy,y)] 计算
eig(subs(fh2,[x,y],[0,0]))
得知正定,所以是极小值点.极小值用
subs(f,[x,y],[0,0])
求得。同理可得(-1/3,-6)为极大值点,其它两个为鞍点。
%Exercise 9(以第一小题为例)
syms y;f=exp(2*y)/(exp(y)+2);
fi=int(f,y)
s=simple(diff(fi)-f)
%Exercise 10
syms x y;f=(x-y)^3*sin(x+2*y);Ix=simple(int(f,y,-x,x)) %Exercise 11(3)
syms x;f=x*log(x^4)*asin(1/x);Ix=int(f,x,1,3);vpa(Ix) %Exercise 12
%1(3)
syms x;solve(5*x^23-6*x^7+8*x^6-5*x^2)
%6
syms a b;s=solve(a-0.7*sin(a)-0.2*cos(b),b-0.7*cos(a)+0.2*sin(b));
s.a,s.b
%Exercise 13
%1(3)
dsolve('D2y-0.01*Dy^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')(解不出)
%1(4)
dsolve('2*D2x-5*Dx+3*x=45*exp(2*t)','x(0)=2','Dx(0)=1','t') %Exercise 14
%6(ii)
ezplot('x^2/4+y^2/9=1')
%6(vi)
ezmesh('2*sin(ph)*cos(th)','2*sin(ph)*sin(th)','2*cos(ph)',[0 pi/2 0 2*pi])
转载请注明出处范文大全网 » matlab数学实验答案
4说明有无穷多解>=1.1).*(x>=1.1).*(x>2)+(t.*t-2.*t+1).*(t>1)+(t.*t-1).*(t>