Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
- 164
- 165
- 166
- 167
- 168
- 169
- 170
- 171
- 172
- 173
- 174
- 175
- 176
- 177
- 178
- 179
- 180
- 181
- 182
- 183
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 202
- 203
- 204
- 205
- 206
- 207
- 208
- 209
- 210
- 211
- 212
- 213
- 214
- 215
- 216
- 217
- 218
- 219
- 220
- 221
- 222
- 223
- 224
- 225
- 226
- 227
- 228
- 229
- 230
- 231
- 232
- 233
- 234
- 235
- 236
- 237
- 238
- 239
- 240
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 259
- 260
- 261
- 262
- 263
- 264
- 265
- 266
- 267
- 268
- 269
- 270
- 271
- 272
- 273
- 274
- 275
- 276
- 277
- 278
- 279
- 280
- 281
- 282
- 283
- 284
- 285
- 286
- 287
- 288
- 289
- 290
- 291
- 292
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 315
- 316
- 317
- 318
- 319
- 320
- 321
- 322
- 323
- 324
- 325
- 326
- 327
- 328
- 329
- 330
- 331
- 332
- 333
- 334
- 335
- 336
- 337
- 338
- 339
- 340
- 341
- 342
- 343
- 344
- 345
- 346
- 347
- 348
- 349
- 350
- 351
- 352
- 353
- 354
- 355
- 356
- 357
- 358
- 359
- 360
- 361
- 362
- 363
- 364
- 365
- 366
- 367
- 368
- 369
- 370
- 371
- 372
- 373
- 374
- 375
- 376
- 377
- 378
- 379
- 380
- 381
- 382
- 383
- 384
- 385
- 386
- 387
- 388
- 389
- 390
- 391
- 392
- 393
- 394
- 395
- 396
- 397
- 398
- 399
- 400
- 401
- 402
- 403
- 404
- 405
- 406
- 407
- 408
- 409
- 410
- 411
- 412
- 413
- 414
- 415
- 416
- 417
- 418
- 419
- 420
- 421
- 422
- 423
- 424
- 425
- 426
- 427
- 428
- 429
- 430
- 431
- 432
- 433
- 434
- 435
- 436
- 437
- 438
- 439
- 440
- 441
- 442
- 443
- 444
- 445
- 446
- 447
- 448
- 449
- 450
- 451
- 452
- 453
- 454
- 455
- 456
- 457
- 458
- 459
- 460
- 461
- 462
- 463
- 464
- 465
- 466
- 467
- 468
- 469
- 470
- 471
- 472
- 473
- 474
- 475
- 476
- 477
- 478
- 479
- 480
- 481
- 482
- 483
- 484
- 485
- 486
- 487
- 488
- 489
- 490
- 491
- 492
- 493
- 494
- 495
- 496
- 497
- 498
- 499
- 500
- 501
- 502
- 503
- 504
- 505
- 506
- 507
- 508
- 509
- 510
- 511
- 512
- 513
- 514
- 515
- 516
- 517
- 518
- 519
- 520
- 521
- 522
- 523
- 524
- 525
- 526
- 527
- 528
- 529
- 530
- 531
- 532
- 533
- 534
- 535
- 536
- 537
- 538
- 539
- 540
- 541
- 542
- 543
- 544
- 545
- 546
- 547
- 548
- 549
- 550
- 551
- 552
- 553
- 554
- 555
- 556
- 557
- 558
- 559
- 560
- 561
- 562
- 563
- 564
- 565
- 566
- 567
- 568
- 569
- 570
- 571
- 572
- 573
- 574
- 575
- 576
- 577
- 578
- 579
- 580
- 581
- 582
- 583
- 584
- 585
- 586
- 587
- 588
- 589
- 590
- 591
- 592
- 593
- 594
- 595
- 596
- 597
- 598
- 599
- 600
- 601
- 602
- 603
- 604
- 605
- 606
- 607
- 608
- 609
- 610
- 611
- 612
- 613
- 614
- 615
- 616
- 617
- 618
- 619
- 620
- 621
- 622
- 623
- 624
- 625
- 626
- 627
- 628
- 629
- 630
- 631
- 632
- 633
- 634
- 635
- 636
- 637
- 638
- 639
- 640
- 641
- 642
- 643
- 644
- 645
- 646
- 647
- 648
- 649
- 650
- 651
- 652
- 653
- 654
- 655
- 656
- 657
- 658
- 659
- 660
- 661
- 662
- 663
- 664
- 665
- 666
- 667
- 668
- 669
- 670
- 671
- 672
- 673
- 674
- 675
- 676
- 677
- 678
- 679
- 680
- 681
- 682
- 683
- 684
- 685
- 686
- 687
- 688
- 689
- 690
- 691
- 692
- 693
- 694
- 695
- 696
- 697
- 698
- 699
- 700
- 701
- 702
- 703
- 704
- 705
- 706
- 707
- 708
- 709
- 710
- 711
- 712
- 713
- 714
- 715
- 716
- 717
- 718
- 719
- 720
- 721
- 722
- 723
- 724
- 725
- 726
- 727
- 728
- 729
- 730
- 731
- 732
- 733
- 734
- 735
- 736
- 737
- 738
- 739
- 740
- 741
- 742
- 743
- 744
- 745
- 746
- 747
- 748
- 749
- 750
- 751
- 752
- 753
- 754
- 755
- 756
- 757
- 758
- 759
- 760
- 761
- 762
- 763
- 764
- 765
- 766
- 767
- 768
- 769
- 770
- 771
- 772
- 773
- 774
- 775
- 776
- 777
- 778
- 779
- 780
- 781
- 782
- 783
- 784
- 785
- 786
- 787
- 788
- 789
- 790
- 791
- 792
- 793
- 794
- 795
- 796
- 797
- 798
- 799
- 800
- 801
- 802
- 803
- 804
- 805
- 806
- 807
- 808
- 809
- 810
- 811
- 812
- 813
- 814
- 815
- 816
- 817
- 818
- 819
- 820
- 821
- 822
- 823
- 824
- 825
- 826
- 827
- 828
- 829
- 830
- 831
- 832
- 833
- 834
- 835
- 836
- 837
- 838
- 839
- 840
- 841
- 842
- 843
- 844
- 845
- 846
- 847
- 848
- 849
- 850
- 851
- 852
- 853
- 854
- 855
- 856
- 857
- 858
- 859
- 860
- 861
- 862
- 863
- 864
- 865
- 866
- 867
- 868
- 869
- 870
- 871
- 872
- 873
- 874
- 875
- 876
- 877
- 878
- 879
- 880
- 881
- 882
- 883
- 884
- 885
- 886
- 887
- 888
- 889
- 890
- 891
- 892
- 893
- 894
- 895
- 896
- 897
- 898
- 899
- 900
- 901
- 902
- 903
- 904
- 905
- 906
- 907
- 908
- 909
- 910
- 911
- 912
- 913
- 914
- 915
- 916
- 917
- 918
- 919
- 920
- 921
- 922
- 923
- 924
- 925
- 926
- 927
- 928
- 929
- 930
- 931
- 932
- 933
- 934
- 935
- 936
- 937
- 938
- 939
- 940
- /***************************************************************/
- /** **/
- /** Tamilyn Tanaka 9373312 **/
- /** Wesley Tsutomu Takai 7990443 **/
- /** **/
- /***************************************************************/
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #include <string.h>
- typedef struct //declaraзгo da estrutura matriz
- {
- int linhas, colunas; //linhas e colunas
- double ** val; // corpo da matriz
- } matriz_t, *matriz;
- matriz matrizNova(int linhas, int colunas);
- void matrizClr(matriz M);
- // exibe na tela a matriz m
- void matrizImp(matriz m);
- //Calcula RKF45, quando f(t,x) й uma funзгo
- void RKF45Funcao(double t0, double x0,double h, double epslon, double tf, int funcao, char* nome,int caso);
- // calcula os valores de f(t, x) do teste 1
- double funcaoTeste1(double t,double x);
- //calcula RKF45, quando F(t,x)=A*X (forma matricial)
- void RKF45Matriz(double t0, matriz A, matriz x0, double h, double epslon, double tf, char* nome, int caso);
- //Caso utilizado para a resoluзгo dos exercicios
- void RKF45MatrizEX(double t0, double x0, double y0, double z0,double h, double epslon, double tf, double c1, double c2, double r, double l,int funcao,char* nome, int caso);
- //multipliaзгo de matrizes
- matriz multiplicaAX(matriz A, matriz X);
- //multiplica um numero constante por uma matriz
- matriz multiplicaCteMatriz(double h,matriz Fx);
- //soma de duas matrizes de mesma dimensгo
- matriz somaMatriz(matriz A, matriz B);
- double maxTl(matriz A);
- double resistorNaoLinear(double V);
- double funcaoVC1(double x0,double x,double y,double z,double R,double C1)
- {
- double resposta;
- resposta= (1/(R*C1))*(y-x) -(1/C1)*resistorNaoLinear(x0);
- return resposta;
- }
- double funcaoVC2(double x,double y,double z,double R,double C2)
- {
- double resposta;
- resposta= (1/(R*C2))*(x-y) +(1/C2)*z;
- return resposta;
- }
- double funcaoI(double x,double y,double z,double R,double L)
- {
- double resposta;
- resposta= -(1/L)*y;
- return resposta;
- }
- int main()
- {
- //declaracao de variaveis
- matriz A,X;
- int escolha, exercicio, teste;
- double t0, x0,tf,epslon, h;
- char* nome="Arquivo";
- //lista de comandos
- printf("escolha:\n 1)exercicios \n 2) testes\n \3)Exercicios extras\n");
- scanf("%d",&escolha);
- while(escolha==1 || escolha==2 || escolha ==3)
- {
- if(escolha==1)
- {
- printf("escolha o numero do exercicio \n" );
- printf("1) r<1500 \n" );
- printf("2) r=1600\n" );
- printf("3) 1500<r<1600\n" );
- printf("4 e 5) R variavel\n" );
- scanf("%d", &exercicio);
- double R=1600, C1=10e-9,C2=100e-9,L=18e-3;
- if(exercicio==1)
- {
- while( R >= 1500 || R <=0)
- {
- printf("Valor de R:" );
- scanf("%lf", &R);
- }
- t0=0;
- h=1e-3;
- epslon=0.001;
- tf=0.05;
- //matrizImp(A);
- nome= "Exercicio1";
- RKF45MatrizEX(t0,-.5,-.2,0,h,epslon,tf,C1,C2,R,L,1,nome,1);
- }
- else if(exercicio==2)
- {
- t0=0;
- h=1e-3;
- epslon=0.001;
- tf=.05;
- //matrizImp(A);
- nome= "Exercicio2";
- RKF45MatrizEX(t0,-.5,-.2,0,h,epslon,tf,C1,C2,R,L,1,nome,1);
- }
- else if(exercicio==3)
- {
- printf("A anбlise do exercicio estб no relatуrio");
- }
- else if(exercicio==4)
- {
- printf("Valor de R:" );
- scanf("%lf", &R);
- while( R <=0)
- {
- printf("Valor de R:" );
- scanf("%lf", &R);
- }
- t0=0;
- h=1e-3;
- epslon=0.001;
- tf=.05;
- //matrizImp(A);
- nome= "Exercicio4";
- RKF45MatrizEX(t0,-.5,-.2,0,h,epslon,tf,C1,C2,R,L,1,nome,1);
- }
- else
- printf("exercicio nao encontrado\n");
- }
- if(escolha==2)
- {
- printf("escolha o numero do teste (de 1 a 3) \n" );
- scanf("%d", &teste);
- if(teste==1)
- {
- //valores fornecidos no enunciado:
- t0=1.05;
- x0=-18.95;
- tf=3;
- epslon=1e-5;
- h=0.1;
- nome ="Teste1";
- RKF45Funcao(t0,x0,h,epslon,tf,1,nome,1);
- }
- else if(teste==2)
- {
- A= matrizNova(4,4);
- A->val[0][0]=A->val[0][3]=A->val[1][1]=A->val[2][1]= A->val[2][2]=A->val[3][3]=-2;
- A->val[0][1]=A->val[0][2]=A->val[1][3]=A->val[2][0]=A->val[2][3]=A->val[3][1]=-1;
- A->val[1][0]=A->val[3][2]=1;
- A->val[1][2]=A->val[3][0]=2;
- X=matrizNova(4,1);
- X->val[0][0]=1;
- X->val[1][0]=1;
- X->val[2][0]=1;
- X->val[3][0]=-1;
- t0=0;
- h=0.1;
- epslon=1e-5;
- tf=2;
- //matrizImp(A);
- nome= "Teste2";
- RKF45Matriz(t0, A, X, h, epslon, tf,nome,2);
- matrizClr(X);
- matrizClr(A);
- }
- else if(teste==3)
- {
- int linhas,i,j;
- printf("Digite o numero de linhas/colunas da matriz quadrada:");
- scanf("%d", & linhas);
- A= matrizNova(linhas, linhas);
- for(i=0; i<linhas; i++)
- {
- for(j=0; j<linhas; j++)
- {
- if(i==j)
- A->val[i][j]=-2;
- else if(i==j+1 || i+1==j)
- A->val[i][j]=1;
- else
- A->val[i][j]=0;
- }
- }
- X=matrizNova(linhas,1);
- for (i=0; i< linhas; i++)
- {
- X->val[i][0]= sin(M_PI*(1+i)/(linhas+1))+sin(linhas*M_PI*(1+i)/(linhas+1));
- }
- matrizImp(X);
- nome ="Teste3";
- t0=0;
- h=0.1;
- epslon=1e-5;
- tf=2;
- RKF45Matriz(t0, A, X, h, epslon, tf,nome,3);
- matrizImp(A);
- matrizClr(A);
- }
- else
- printf("teste nao encontrado \n");
- }
- if(escolha==3)
- {
- printf("Escolha o exercicio extra 1 ou 2\n" );
- scanf("%d", &exercicio);
- char* nome="ExercicioExtra1";
- double v10=-.5, v20=-.2, i0=0;
- double R=20000, C1=10e-9,C2=100e-9,L=18e-3;
- if(exercicio ==1 )
- {
- char* nome="ExercicioExtra1";
- printf("Valor de vc10:" );
- scanf("%lf", &v10);
- printf("Valor de vc20:" );
- scanf("%lf", &v20);
- printf("Valor de i0:" );
- scanf("%lf", &i0);
- printf("Valor de R:" );
- scanf("%lf", &R);
- while( R <=0)
- {
- printf("Valor de vc10:" );
- scanf("%lf", &v10);
- printf("Valor de vc20:" );
- scanf("%lf", &v20);
- printf("Valor de i0:" );
- scanf("%lf", &i0);
- printf("Valor de R:" );
- scanf("%lf", &R);
- }
- t0=0;
- h=1e-3;
- epslon=1e-3;
- tf=0.05;
- //matrizImp(A);
- nome= "Exercicio1";
- RKF45MatrizEX(t0,v10,v20,i0,h,epslon,tf,C1,C2,R,L,exercicio,nome,1);
- }
- else if(exercicio ==2 )
- {
- char* nome="ExercicioExtra2";
- printf("Valor de L:" );
- scanf("%lf", &L);
- printf("Valor de C1:" );
- scanf("%lf", &C1);
- printf("Valor de C2:" );
- scanf("%lf", &C2);
- printf("Valor de R:" );
- scanf("%lf", &R);
- t0=0;
- h=1e-3;
- epslon=1e-3;
- tf=0.05;
- //matrizImp(A);
- nome= "Exercicio1";
- RKF45MatrizEX(t0,v10,v20,i0,h,epslon,tf,C1,C2,R,L,1,nome,1);
- }
- else
- {
- printf("nao encontrado\n");
- }
- }
- printf("escolha:\n 1)exercicios \n 2) testes\n 3)Exercicios extras\n");
- scanf("%d",&escolha);
- }
- printf("Fim do programa");
- return 0;
- }
- matriz matrizNova(int linhas, int colunas)
- {
- int i;
- matriz nova = malloc(sizeof(matriz_t));
- nova->val = malloc(sizeof(double) * linhas);
- nova->val[0] = calloc(sizeof(double), linhas * colunas);
- for (i = 0; i < linhas; i++)
- {
- nova->val[i] = nova->val[0] + colunas * i;
- }
- nova->linhas = linhas;
- nova->colunas = colunas;
- return nova;
- }
- void matrizClr(matriz M)
- {
- free(M->val[0]);
- free(M->val);
- free(M);
- }
- void matrizImp(matriz m) // funзгo para imprimir a matriz por varredura
- {
- int i,j;
- for(i = 0; i < m->linhas; i++)
- {
- for (j = 0; j < m->colunas; j++)
- {
- printf(" %e", m->val[i][j]);
- }
- printf("\n");
- }
- printf("\n");
- }
- void RKF45Funcao(double t0, double x0,double h, double epslon, double tf, int funcao,char* nome, int caso)
- {
- double k1,k2,k3,k4,k5,k6;
- double xQuartaOrdem, xQuintaOrdem;
- double t, x, tl1=1,h2;
- double c=2;
- t=t0;
- x=x0;
- char* prefixo_do_nome =nome;
- FILE *arq;
- char fname[200];
- /* monta o nome do arquivo */
- strcpy(fname, prefixo_do_nome);
- strcat(fname,".csv");
- /*abre o o arquivo para gravacao */
- arq = fopen(fname, "w");
- if (arq == NULL)
- {
- printf("Erro na abertura do arquivo %s\n\n",fname);
- return 0;
- }
- //calculo de k1,k2,k3,k4,k5 e k6 de acordo com a funзгo de cada exercicio
- fprintf(arq,"t x h Erro\n");
- if(funcao==1)
- {
- while (t0 +h<=tf)
- {
- while (fabs(tl1) > epslon) //ajuste do passo e calculo dos ks
- {
- //utiliza a funcao fornecida no teste 1
- k1=h*funcaoTeste1(t0,x0);
- //printf("k1= %f \n", k1);
- t=t0+h*1/4;
- x=x0+k1*1/4;
- k2=h*funcaoTeste1(t,x);
- //printf("k2= %f \n", k2);
- t=t0+h*3/8;
- x=x0+k1*3/32+k2*9/32;
- k3=h*funcaoTeste1(t,x);
- //printf("k3= %f \n", k3);
- t=t0+h*12/13;
- x=x0+k1*1932/2197-k2*7200/2197+k3*7296/2197;
- k4=h*funcaoTeste1(t,x);
- //printf("k4= %f \n", k4);
- t=t0+h;
- x=x0+k1*439/216-8*k2+k3*3680/513-k4*845/4104;
- k5=h*funcaoTeste1(t,x);
- //printf("k5= %f \n", k5);
- t=t0+h*1/2;
- x=x0-k1*8/27+2*k2-k3*3544/2565+k4*1859/4104-k5*11/40;
- k6=h*funcaoTeste1(t,x);
- //printf("k6= %f \n", k6);
- // calculo de Xi+1 de quarta ordem e de X~i+1 de quinta ordem com os valores de k1,k2,k3,k4,k5 e k6 obtidos anteriormnete
- xQuartaOrdem=x0+k1*25/216+k3*1408/2565+k4*2197/4104-k5*1/5;
- //printf("x quarta ordem= %f \n", xQuartaOrdem);
- xQuintaOrdem=x0+k1*16/135+k3*6656/12825+k4*28561/56430-k5*9/50+k6*2/55;
- //printf("x quinta ordem= %f \n",xQuintaOrdem);
- tl1= (xQuintaOrdem-xQuartaOrdem)/h;
- h2=h;
- h =h*pow((epslon*h)/(c*fabs(tl1*h)),0.25);
- if (h> 1)
- {
- h=1.0;
- }
- if (h < 0.001)
- {
- h=0.001;
- }
- //printf("tl1= %e \nh2= %e \n\n", tl1,h);
- }
- x0= xQuartaOrdem;
- t0=t0+h2;
- t=t0;
- tl1=1;
- fprintf(arq,"%e %e %e %e\n",t0, x0,h,(t +(1/(1-t))-x0));
- printf("%e %e %e %e\n",t0, x0,h,(t +(1/(1-t))-x0));
- if(t0+h>tf && t0!=tf)
- {
- h=tf-t0;
- }
- }
- }
- else
- {
- printf("erro ao tentar calcular o metodo RKF45");
- }
- fclose(arq); /* fecha arquivo */
- }
- double funcaoTeste1(double t,double x)
- {
- double resposta;
- resposta= 1+((x-t)*(x-t));
- return resposta;
- }
- void RKF45Matriz(double t0, matriz A, matriz x0, double h, double epslon, double tf,char* nome,int caso)
- {
- int i;
- matriz k1,k2,k3,k4,k5,k6,X,Fx, temporario;
- matriz xQuartaOrdem, xQuintaOrdem,tl1;
- double t,x, max=1, h2;
- double c=2;
- k1=matrizNova(x0->linhas,x0->colunas);
- k2=matrizNova(x0->linhas,x0->colunas);
- k3=matrizNova(x0->linhas,x0->colunas);
- k4=matrizNova(x0->linhas,x0->colunas);
- k5=matrizNova(x0->linhas,x0->colunas);
- k6=matrizNova(x0->linhas,x0->colunas);
- X=matrizNova(x0->linhas,x0->colunas);
- Fx=matrizNova(x0->linhas,x0->colunas);
- temporario=matrizNova(x0->linhas,x0->colunas);
- xQuartaOrdem=matrizNova(x0->linhas,x0->colunas);
- xQuintaOrdem=matrizNova(x0->linhas,x0->colunas);
- char* prefixo_do_nome =nome;
- FILE *arq;
- char fname[200];
- /* monta o nome do arquivo */
- strcpy(fname, prefixo_do_nome);
- strcat(fname,".csv");
- /*abre o o arquivo para gravacao */
- arq = fopen(fname, "w");
- if (arq == NULL)
- {
- printf("Erro na abertura do arquivo %s\n\n",fname);
- return 0;
- }
- if (caso==2)
- fprintf(arq,"x1 x2 x3 x4 t h Errox1 Errox2 Errox3 Errox4\n");
- if (caso==3)
- {
- for(int j=1; j <= X->linhas; j++)
- {
- fprintf(arq,"x%d ",j);
- }
- fprintf(arq,"t h ");
- for(int j=1; j <= X->linhas; j++)
- {
- fprintf(arq,"Errox%d ",j);
- }
- fprintf(arq,"\n");
- }
- while (t0 +h <= tf)
- {
- /* code */
- while (fabs(max) > epslon)
- {
- // calculo da matriz k1
- for(i=0; i<X->linhas; i++)
- {
- X->val[i][0]=x0->val[i][0];
- }
- Fx=multiplicaAX(A,X);
- k1=multiplicaCteMatriz(h,Fx);
- //printf("k1= \n");
- //matrizImp(k1);
- // calculo da matriz k2
- double a =1;
- temporario=multiplicaCteMatriz(a/4,k1);
- X=somaMatriz(x0,temporario);
- Fx=multiplicaAX(A,X);
- k2=multiplicaCteMatriz(h,Fx);
- //printf("k2= \n");
- //matrizImp(k2);
- // calculo da matriz k3
- temporario=multiplicaCteMatriz(a*3/32,k1);
- X=somaMatriz(x0,temporario);
- temporario=multiplicaCteMatriz(a*9/32,k2);
- X=somaMatriz(X,temporario);
- Fx=multiplicaAX(A,X);
- k3=multiplicaCteMatriz(h,Fx);
- //printf("k3= \n");
- //matrizImp(k3);
- a =1;
- // calculo da matriz k4
- temporario=multiplicaCteMatriz(a*1932/2197,k1);
- X=somaMatriz(x0,temporario);
- temporario=multiplicaCteMatriz(-a*7200/2197,k2);
- X=somaMatriz(X,temporario);
- temporario=multiplicaCteMatriz(a*7296/2197,k3);
- X=somaMatriz(X,temporario);
- Fx=multiplicaAX(A,X);
- k4=multiplicaCteMatriz(h,Fx);
- // printf("k4= \n");
- // matrizImp(k4);
- // calculo da matriz k5
- temporario=multiplicaCteMatriz(a*439/216,k1);
- X=somaMatriz(x0,temporario);
- temporario=multiplicaCteMatriz(-a*8,k2);
- X=somaMatriz(X,temporario);
- temporario=multiplicaCteMatriz(a*3680/513,k3);
- X=somaMatriz(X,temporario);
- temporario=multiplicaCteMatriz(-a*845/4104,k4);
- X=somaMatriz(X,temporario);
- Fx=multiplicaAX(A,X);
- k5=multiplicaCteMatriz(h,Fx);
- // printf("k5= \n");
- // matrizImp(k5);
- // calculo da matriz k6
- temporario=multiplicaCteMatriz(-a*8/27,k1);
- X=somaMatriz(x0,temporario);
- temporario=multiplicaCteMatriz(a*2,k2);
- X=somaMatriz(X,temporario);
- temporario=multiplicaCteMatriz(-a*3544/2565,k3);
- X=somaMatriz(X,temporario);
- temporario=multiplicaCteMatriz(a*1859/4104,k4);
- X=somaMatriz(X,temporario);
- temporario=multiplicaCteMatriz(-a*11/40,k5);
- X=somaMatriz(X,temporario);
- Fx=multiplicaAX(A,X);
- k6=multiplicaCteMatriz(h,Fx);
- // printf("k6= \n");
- // matrizImp(k6);
- //calcula x de quarta ordem
- temporario=multiplicaCteMatriz(a*25/216,k1);
- xQuartaOrdem=somaMatriz(x0,temporario);
- temporario=multiplicaCteMatriz(a*1408/2565,k3);
- xQuartaOrdem=somaMatriz(xQuartaOrdem,temporario);
- temporario=multiplicaCteMatriz(a*2197/4104,k4);
- xQuartaOrdem=somaMatriz(xQuartaOrdem,temporario);
- temporario=multiplicaCteMatriz(-a*0.2,k5);
- xQuartaOrdem=somaMatriz(xQuartaOrdem,temporario);
- // printf("x quarta ordem=\n");
- // matrizImp(xQuartaOrdem);
- //calcula x de quinta ordem
- temporario=multiplicaCteMatriz(a*16/135,k1);
- xQuintaOrdem=somaMatriz(x0,temporario);
- temporario=multiplicaCteMatriz(a*6656/12825,k3);
- xQuintaOrdem=somaMatriz(xQuintaOrdem,temporario);
- temporario=multiplicaCteMatriz(a*28561/56430,k4);
- xQuintaOrdem=somaMatriz(xQuintaOrdem,temporario);
- temporario=multiplicaCteMatriz(-a*9/50,k5);
- xQuintaOrdem=somaMatriz(xQuintaOrdem,temporario);
- temporario=multiplicaCteMatriz(a*2/55,k6);
- xQuintaOrdem=somaMatriz(xQuintaOrdem,temporario);
- // printf("x quinta ordem=\n");
- // matrizImp(xQuintaOrdem);
- max= (maxTl(somaMatriz(xQuintaOrdem, multiplicaCteMatriz(-1,xQuartaOrdem))))/h;
- // printf("tal=\n");
- // matrizImp(somaMatriz(xQuintaOrdem, multiplicaCteMatriz(-1,xQuartaOrdem)));
- h2=h;
- h =h*pow((epslon*h)/(c*fabs(max*h)),0.25);
- if (h> 1)
- {
- h=1.0;
- }
- if (h < 0.001)
- {
- h=0.001;
- }
- printf("h: %e t: %e max|X -X'|: %e\n",h,t0,fabs(max*h));
- }
- x0= xQuartaOrdem;
- t0=t0+h2;
- t=t0;
- max=1;
- int i,j;
- for(i = 0; i < x0->colunas; i++)
- {
- for (j = 0; j < x0->linhas; j++)
- {
- fprintf(arq,"%e ", x0->val[i][j]);
- }
- fprintf(arq,"%e %e ",t0,h);
- }
- //printf("%e %e %e\n",t0,h,exp(-t0)*sin( t0) + exp(-3*t0)*cos(3*t0)-x0->val[0][0]);
- if (caso==2)
- {
- fprintf(arq,"%e %e %e %e\n",exp(-t0)*sin( t0) + exp(-3*t0)*cos(3*t0)-x0->val[0][0],exp(-t0)*cos( t0) + exp(-3*t0)*sin(3*t0)-x0->val[1][0],-exp(-t0)*sin( t0) + exp(-3*t0)*cos(3*t0)-x0->val[2][0],-exp(-t0)*cos( t0) + exp(-3*t0)*sin(3*t0)-x0->val[3][0]);
- }
- if (caso==3)
- {
- for(int i=0; i < x0->linhas; i++)
- {
- double del1 =(-2)*(1-cos(M_PI/(x0->linhas+1)))*t0;
- double del2 =(-2)*(1-cos((M_PI*x0->linhas)/(x0->linhas+1)))*t0;
- double e = exp(del1)*sin(M_PI*(1+i)/(x0->linhas+1)) +exp(del2)*sin(x0->linhas*M_PI*(1+i)/(x0->linhas+1)) -x0->val[i][0];
- fprintf(arq,"%e ",e);
- }
- fprintf(arq,"%\n");
- //fprintf(arq,"%e %e %e %e\n",exp(-t0)*sin( t0) + exp(-3*t0)*cos(3*t0)-x0->val[0][0],exp(-t0)*cos( t0) + exp(-3*t0)*sin(3*t0)-x0->val[1][0],-exp(-t0)*sin( t0) + exp(-3*t0)*cos(3*t0)-x0->val[2][0],-exp(-t0)*cos( t0) + exp(-3*t0)*sin(3*t0)-x0->val[3][0]);
- }
- if(t0+h>tf && t0!=tf)
- {
- h=tf-t0;
- }
- }
- fclose(arq);
- matrizClr(k1);
- matrizClr(k2);
- matrizClr(k3);
- matrizClr(k4);
- matrizClr(k5);
- matrizClr(k6);
- matrizClr(X);
- matrizClr(Fx);
- matrizClr(temporario);
- matrizClr(xQuartaOrdem);
- matrizClr(xQuintaOrdem);
- //matrizClr(tl1);
- }
- void RKF45MatrizEX(double t0, double x0, double y0, double z0,double h, double epslon, double tf, double c1, double c2, double r, double l,int funcao,char* nome, int caso)
- {
- double xk1,xk2,xk3,xk4,xk5,xk6;
- double xQuartaOrdem, xQuintaOrdem;
- double yk1,yk2,yk3,yk4,yk5,yk6;
- double yQuartaOrdem, yQuintaOrdem;
- double zk1,zk2,zk3,zk4,zk5,zk6;
- double zQuartaOrdem, zQuintaOrdem;
- double t=t0, x=x0,y=y0,z=z0, tl1=1,h2;
- double c=2;
- int i=0;
- char* prefixo_do_nome =nome;
- FILE *arq;
- char fname[200];
- /* monta o nome do arquivo */
- strcpy(fname, prefixo_do_nome);
- strcat(fname,".csv");
- /*abre o o arquivo para gravacao */
- arq = fopen(fname, "w");
- if (arq == NULL)
- {
- printf("Erro na abertura do arquivo %s\n\n",fname);
- return 0;
- }
- //calculo de k1,k2,k3,k4,k5 e k6 de acordo com a funзгo de cada exercicio
- fprintf(arq,"t Vc1 Vc2 IL h\n");
- fprintf(arq,"%e %e %e %e %e\n",t0, x0,y0,z0,h);
- if(funcao==1)
- {
- while (t0 +h<=tf)
- {
- i++;
- while (fabs(tl1) > epslon) //ajuste do passo e calculo dos ks
- {
- //utiliza a funcao fornecida no teste 1
- xk1=h*funcaoVC1(x0,x0,y0,z0,r,c1);
- yk1=h*funcaoVC2(x0,y0,z0,r,c2);
- zk1=h*funcaoI(x0,y0,z0,r,l);
- // printf("xk1= %e \n", xk1);
- // printf("yk1= %e \n", yk1);
- // printf("zk1= %e \n\n", zk1);
- t=t0+h*1/4;
- x=x0+xk1*1/4;
- y=y0+yk1*1/4;
- z=z0+zk1*1/4;
- // printf("x= %e \n", x);
- // printf("y= %e \n", y);
- // printf("z= %e \n\n", z);
- xk2=h*funcaoVC1(x0,x,y,z,r,c1);
- yk2=h*funcaoVC2(x,y,z,r,c2);
- zk2=h*funcaoI(x,y,z,r,l);
- // printf("xk2= %e \n", xk2);
- // printf("yk2= %e \n", yk2);
- // printf("zk2= %e \n\n", zk2);
- t=t0+h*3/8;
- x=x0+xk1*3/32+xk2*9/32;
- y=y0+yk1*3/32+yk2*9/32;
- z=z0+zk1*3/32+zk2*9/32;
- xk3=h*funcaoVC1(x0,x,y,z,r,c1);
- yk3=h*funcaoVC2(x,y,z,r,c2);
- zk3=h*funcaoI(x,y,z,r,l);
- // printf("xk3= %e \n", xk3);
- // printf("yk3= %e \n", yk3);
- // printf("zk3= %e \n\n", zk3);
- t=t0+h*12/13;
- x=x0+xk1*1932/2197-xk2*7200/2197+xk3*7296/2197;
- y=y0+yk1*1932/2197-yk2*7200/2197+yk3*7296/2197;
- z=z0+zk1*1932/2197-zk2*7200/2197+zk3*7296/2197;
- xk4=h*funcaoVC1(x0,x,y,z,r,c1);
- yk4=h*funcaoVC2(x,y,z,r,c2);
- zk4=h*funcaoI(x,y,z,r,l);
- // printf("xk4= %e \n", xk4);
- // printf("yk4= %e \n", yk4);
- // printf("zk4= %e \n\n", zk4);
- t=t0+h;
- x=x0+xk1*439/216-8*xk2+xk3*3680/513-xk4*845/4104;
- y=y0+yk1*439/216-8*yk2+yk3*3680/513-yk4*845/4104;
- z=z0+zk1*439/216-8*zk2+zk3*3680/513-zk4*845/4104;
- xk5=h*funcaoVC1(x0,x,y,z,r,c1);
- yk5=h*funcaoVC2(x,y,z,r,c2);
- zk5=h*funcaoI(x,y,z,r,l);
- // printf("xk5= %e \n", xk5);
- // printf("yk5= %e \n", yk5);
- // printf("zk5= %e \n\n", zk5);
- t=t0+h*1/2;
- x=x0-xk1*8/27+2*xk2-xk3*3544/2565+xk4*1859/4104-xk5*11/40;
- y=y0-yk1*8/27+2*yk2-yk3*3544/2565+yk4*1859/4104-yk5*11/40;
- z=z0-zk1*8/27+2*zk2-zk3*3544/2565+zk4*1859/4104-zk5*11/40;
- xk6=h*funcaoVC1(x0,x,y,z,r,c1);
- yk6=h*funcaoVC2(x,y,z,r,c2);
- zk6=h*funcaoI(x,y,z,r,l);
- // printf("xk6= %e \n", xk6);
- // printf("yk6= %e \n", yk6);
- // printf("zk6= %e \n\n", zk6);
- // calculo de Xi+1 de quarta ordem e de X~i+1 de quinta ordem com os valores de k1,k2,k3,k4,k5 e k6 obtidos anteriormnete
- xQuartaOrdem=x0+xk1*25/216+xk3*1408/2565+xk4*2197/4104-xk5*1/5;
- yQuartaOrdem=y0+yk1*25/216+yk3*1408/2565+yk4*2197/4104-yk5*1/5;
- zQuartaOrdem=z0+zk1*25/216+zk3*1408/2565+zk4*2197/4104-zk5*1/5;
- // printf("x quarta ordem= %e \n", xQuartaOrdem);
- // printf("y quarta ordem= %e \n", yQuartaOrdem);
- // printf("z quarta ordem= %e \n\n", zQuartaOrdem);
- xQuintaOrdem=x0+xk1*16/135+xk3*6656/12825+xk4*28561/56430-xk5*9/50+xk6*2/55;
- yQuintaOrdem=y0+yk1*16/135+yk3*6656/12825+yk4*28561/56430-yk5*9/50+yk6*2/55;
- zQuintaOrdem=z0+zk1*16/135+zk3*6656/12825+zk4*28561/56430-zk5*9/50+zk6*2/55;
- // printf("x quinta ordem= %e \n", xQuintaOrdem);
- // printf("y quinta ordem= %e \n", yQuintaOrdem);
- // printf("z quinta ordem= %e \n\n", zQuintaOrdem);
- if(fabs(xQuintaOrdem-xQuartaOrdem) > fabs(yQuintaOrdem-yQuartaOrdem) && fabs(xQuintaOrdem-xQuartaOrdem) > fabs(zQuintaOrdem-zQuartaOrdem) )
- {
- tl1= (xQuintaOrdem-xQuartaOrdem)/h;
- }
- else if(fabs(yQuintaOrdem-yQuartaOrdem) > fabs(xQuintaOrdem-xQuartaOrdem) && fabs(yQuintaOrdem-yQuartaOrdem) > fabs(zQuintaOrdem-zQuartaOrdem) )
- {
- tl1= (yQuintaOrdem-yQuartaOrdem)/h;
- }
- else
- {
- tl1= (zQuintaOrdem-zQuartaOrdem)/h;
- }
- h2=h;
- h =h*pow((epslon*h)/(c*fabs(tl1*h)),0.25);
- if (h> 1)
- {
- h=1.0;
- }
- if (h < 1e-10)
- {
- h=1e-8;
- }
- // printf("tl1= %e \nh2= %e \n\n", tl1,h);
- }
- x0= xQuartaOrdem;
- y0= yQuartaOrdem;
- z0= zQuartaOrdem;
- t0=t0+h2;
- t=t0;
- tl1=1;
- fprintf(arq,"%e %e %e %e %e\n",t0, x0,y0,z0,h);
- printf("%e %e %e %e %e\n",t0, x0,y0,z0,h);
- if(t0+h>tf && t0!=tf)
- {
- h=tf-t0;
- }
- }
- }
- else
- {
- printf("erro ao tentar calcular o metodo RKF45");
- }
- fclose(arq); /* fecha arquivo */
- }
- double maxTl(matriz A)
- {
- double max= A->val[0][0];
- int i;
- for (i=0; i < A->linhas; i++)
- {
- if(fabs(A->val[i][0]) > fabs(max) )
- {
- max = A->val[i][0];
- }
- }
- return max;
- }
- matriz multiplicaAX(matriz A, matriz X)
- {
- int i, j, k;
- matriz resposta;
- resposta=matrizNova(A->linhas,X->colunas);
- //zerando a matriz resposta
- for(i=0; i< A->linhas; i++)
- {
- for(j=0; j< X->colunas; j++)
- {
- resposta->val[i][j]=0;
- }
- }
- if(A->colunas==X->linhas)
- {
- //multiplicaзгo de matrizes
- for (i=0; i<A->linhas; i++)
- {
- for (j=0; j<X->colunas; j++)
- {
- for (k=0; k<A->colunas; k++)
- {
- resposta->val[i][j] = resposta->val[i][j] + (A->val[i][k] * X->val[k][j]);
- }
- }
- }
- }
- else
- printf("erro ao tentar efetuar o calculo A*X \n");
- return resposta;
- }
- matriz multiplicaCteMatriz(double h,matriz Fx)
- {
- int i, j;
- matriz resposta;
- resposta=matrizNova(Fx->linhas, Fx->colunas);
- for(i=0; i<Fx->linhas; i++)
- {
- for(j=0; j<Fx->colunas; j++)
- {
- resposta->val[i][j]=h*Fx->val[i][j];
- }
- }
- return resposta;
- }
- matriz somaMatriz(matriz A, matriz B)
- {
- int i,j;
- matriz resposta;
- resposta=matrizNova(A->linhas, A->colunas);
- if(A->linhas==B->linhas && A->colunas==B->colunas)
- {
- for(i=0; i<A->linhas; i++)
- {
- for(j=0; j<A->colunas; j++)
- {
- resposta->val[i][j]=A->val[i][j]+B->val[i][j];
- }
- }
- }
- else
- printf("erro ao tentar somar as matrizes");
- return resposta;
- }
- double resistorNaoLinear(double V)
- {
- double E, Emax, Ga, Gb, Gc, Gv,a=1;
- E=1.17391304;
- Emax=8.1818;
- Ga= -(1e-3)*50/66 ;
- Gb= -(1e-3)*9/22;
- Gc= 4.591*1e-3;
- if(V<=(-Emax)) // < -8.18
- Gv=Gc*V + Emax*(Gc-Gb)+E*(Gb-Ga);
- else if(V<=(-E)) // -8.18 < V <= 1.17
- Gv=Gb*V+(Gb-Ga)*E;
- else if( V<E) //-1.17 < V < 1.17
- Gv=Ga*V;
- else if(V<Emax)
- Gv=Gb*V+(Ga-Gb)*E;
- else
- {
- Gv=Gc*V+Emax*(Gb-Gc)+E*(Ga-Gb);
- }
- return Gv;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement