123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156 |
- #define _CRT_SECURE_NO_WARNINGS
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #include <Windows.h>
- #include "header.h"
- vector* find_initials(double A, double alpha_0, double alpha_1) {
- vector u = { alpha_0*A/(pow(alpha_0,2)+ pow(alpha_1,2)),alpha_1*A / (pow(alpha_0,2) + pow(alpha_1,2)) };
- vector v = { alpha_1,-alpha_0 };
- vector* ans= malloc(sizeof(vector)*2);
- ans[0] = u;
- ans[1] = v;
- return ans;
- }
- vector* euler_method(vector y0, double a, double b, vector(*F)(double, vector), int n) {
- //(x[0], y0) - ïåðâàÿ, èìåþùàÿñÿ ó íàñ òî÷êà òàáëè÷íîé ôóíêöèè
- double h = (b - a) / (n - 1);
- vector* y = malloc(sizeof(vector) * n);
- if (y == NULL) {
- exit(16);
- }
- y[0] = y0;
-
- for (int i = 0;i < n - 1;i++) {
- vector res = F(a + i * h, y[i]);
- y[i + 1].y1 = y[i].y1 + h * res.y1;
- y[i + 1].y2 = y[i].y2 + h * res.y2;
- }
- return y;
- }
- double* solve(vector* initials,double betta_0,double betta_1, double B,double a,double b,vector(*F)(double, vector), vector(*F_0)(double, vector),int n) {
- vector u0 = initials[0];
- vector v0 = initials[1];
- vector* u = euler_method(u0, a, b, F, n);
- vector* v = euler_method(v0, a, b, F_0, n);
- double* y = malloc(sizeof(double) * n);
- double C = (B - betta_0 * u[n - 1].y1 - betta_1 * u[n - 1].y2) / (betta_0 * v[n - 1].y1 + betta_1 * v[n - 1].y2);
- for (int i = 0;i < n;i++) {
-
- y[i] = u[i].y1 + C * v[i].y1;
- }
- return y;
- }
- double** solve_runge(vector* initials, double betta_0, double betta_1, double B, double a, double b, vector(*F)(double, vector), vector(*F_0)(double, vector), double eps, int* n_, int* len_h){
-
- printf("eps = %e\n", eps);
- double* x = malloc(sizeof(double));
- long double* h_all = malloc(sizeof(long double));// ðåãèñòðèðóåò âñå çíà÷åíèÿ øàãîâ
- long double* h_made = malloc(sizeof(long double));// ðåãèñòðèðóåò çíà÷åíèÿ øàãà, ñîîòâåòñòâóþùèå ñäåëàíûì
- long double h = (b - a) / 2;
- x[0] = a;
- h_all[0] = h;
- vector res;
- double c1,c2;
- vector *u1 = malloc(sizeof(vector)*2), * v1 = malloc(sizeof(vector) * 2), u2, v2;
- u1[0] = initials[0];
- v1[0] = initials[1];
- u2 = initials[0];
- v2 = initials[1];
- int i = 1, k = 1, n = 10000;
- vector* u = malloc(sizeof(vector));
- vector* v = malloc(sizeof(vector));
- u[0] = initials[0];
- v[0] = initials[1];
- eps = pow(eps, 3.0/2);
- double difu,difv,dif;
- while (x[i - 1] + h <= b) {
- h_all = realloc(h_all, sizeof(long double) * (k + 1));
- h_all[k] = h;
- u1[1] = euler_method(u1[0], x[i - 1], x[i - 1] + h, F, 3)[2];
- v1[1] = euler_method(v1[0], x[i - 1], x[i - 1] + h, F_0, 3)[2];
- u2 = euler_method(u1[0], x[i - 1], x[i - 1] + h, F, 5)[4];
- v2 = euler_method(v1[0], x[i - 1], x[i - 1] + h, F_0, 5)[4];
- difu = fabs(u1[1].y1 - u2.y1);
- difv = fabs(v1[1].y1 - v2.y1);
-
- //printf("difv = %e, v1[1].y1 = %e, v2.y1 = %e, i = %i,h = %e\n", difv, v1[1].y1, v2.y1,i,h);
-
- if (difu > eps/100 || difv > eps/100) {
- h /= 2;
- }
- else if (difu < eps / 500 && difv < eps / 500) {
- x = realloc(x, sizeof(double) * (i + 1));
- x[i] = x[i - 1] + h;
- u = realloc(u, sizeof(vector) * (i + 1));
- u[i] = u2;
- v = realloc(v, sizeof(vector) * (i + 1));
- v[i] = v2;
- h_made = realloc(h_made, sizeof(long double) * (i + 1));
- h_made[i - 1] = h;
- i += 1;
- u1[0] = u2;
- v1[0] = v2;
- if (x[i] + h * 2 <= b) {
- h *= 2;
- }
- }
- else {
- x = realloc(x, sizeof(double) * (i + 1));
- x[i] = x[i - 1] + h;
- u = realloc(u, sizeof(vector) * (i + 1));
- u[i] = u2;
- v = realloc(v, sizeof(vector) * (i + 1));
- v[i] = v2;
- h_made = realloc(h_made, sizeof(long double) * (i + 1));
- h_made[i-1] = h;
- i += 1;
- u1[0] = u2;
- v1[0] = v2;
- }
- k++;
- }
- h_all[k - 1] = h;
- h_made[i - 1] = h;
- *len_h = k;
- *n_ = i;
- double C = (B - betta_0 * u[i - 1].y1 - betta_1 * u[i - 1].y2) / (betta_0 * v[i - 1].y1 + betta_1 * v[i - 1].y2);
- double* y = malloc(sizeof(double) * i);
- for (int l = 0;l < i;l++) {
- y[l] = u[l].y1 + C * v[l].y1;
- }
- double** ans = malloc(sizeof(double*) * 4);
- ans[0] = x;
- ans[1] = y;
- ans[2] = h_all;
- ans[3] = h_made;
- return ans;
- }
|