Edp
Edp
Edp
Introducci
on
Existen diferentes metodo para resolver Ecuaciones Diferenciales Parciales(EDP) y Ecuaciones Diferenciales Ordinarias(EDO), entre ellos se encuentra el Metodo de las Diferencias Finitas(MDF), el
cual consiste en introducir una malla 1 sobre una region y aproximar las derivadas del problema
planteado por medio de tecnicas de aproximacion de las derivadas como lo es, por ejemplo, la descomposicion en serie de Taylor de la funcion. Ahora tomemos la descomposicion en serie de Taylor
de la funcion U en los puntos x + h y x h, despues sumamos como es habitual en estos casos para
hallar una expresion tanto de la primera como de la segunda derivada de la funcion, entre otras, as
tenemos lo siguiente
U (x + h) U (x) + hU (x) +
h 2
h 2
U (x) +
U (x) + . . .
2!
3!
h 2
h 2
U (x)
U (x) + . . .
2!
3!
U (x + h) + U (x h) = 2U (x) + h2 U (x) + O(h4 ),
U (x h) U (x) hU (x) +
(1)
(2)
(3)
Sea U (x, t) una funcion de x y de t, subdividiremos el plano x t en un conjunto de rectangulos iguales de lados
x y t respectivamente, de tal manera que cada punto de los ejes OX y OY est
an definidos respectivamente por
xi = ih, h = x, i = . . . 1, 0, 1, . . . , yj = jk, k = y, j = 0, 1, . . . , as cada punto de la grilla tendr
a coordenadas
(xi , yj ) = (ih, jk, ).Ver Fig.1.
1
[U (x + h) 2U (x) + U (x h)]
h2
(4)
con un error de h2 .
Ahora de (2) restemos (3), obtenemos la denominada diferencia central
U (x) =
1
[U (x + h) U (x h)]
2h
(5)
1
[U (x + h) U (x)]
h
(6)
1
[U (x) U (x h)]
h
(7)
A continuacion veamos la aplicacion del MDF en una EDO con un problema de Cauchy:
(
u (x) + u(x)
= f (x),
u(a) = , u(b) =
a x b,
= cte
(8)
h = xn+1 xn
(9)
(10)
esto es un sistema de N 1 ecuaciones algebraicas cuyas incognitas son las soluciones aproximadas
en cada retculo de la malla, el n
umero de soluciones un , 0 n N es igual a N + 1, es decir, una
mas que el n
umero de ecuaciones en el sistema (8), faltaran dos ecuaciones que pueden ser suplidas
por las condiciones en los extremos
u0 = u(x0 ) = a,
uN = u(xN ) = b
(11)
(12)
Con el sistema de ecuaciones completo, podemos resolverlo por diferentes metodos numericos(Gausssin Pivoteo, por ejemplo).
2
Formulaci
on del problema
Un ejemplo clasico de una EDP del tipo parabolico es la ecuacion de la distribucion de calor(Difusion)
en una barra de longitud L con las particularidades siguientes:
La barra esta aislada totalmente en sus extremos.
La distribucion de la temperatura es uniforme en las secciones rectas.
El material que compone a la barra es homogeneo.
La ecuacion de distribucion de calor tiene la siguiente forma
U
2U
= 2 , 0 < x < L, t > 0, cte.
t
x
(13)
donde la solucion a esta ecuacion es la temperatura U a una distancia X del extremo de la barra
delgada despues de un tiempo T , la solucion es una funcion que depende de dos variables U (x, t),
es el coeficiente de distribucion de calor.
El M
etodo Explicito
k
(ui+1,j 2ui,j + ui1,j ) + ui,j
h2
(15)
k
, entonces (15) toma la forma
h2
ui,j+1 = rui1,j + (1 2r)ui,j + rui+1,j
(16)
ui,j+1
ui1,j
ui,j
ui+1,j
i 1, j + 1
i, j + 1
i + 1, j + 1
i 1, j
i, j
i + 1, j
(17)
(18)
(19)
donde (18) son las condiciones en la frontera, (19) son las denominadas condiciones iniciales.
La solucion analtica para el problema descrito anteriormente es
n
8 X 1
(sin
)(sin nx)exp(n2 2 t)
U (x, t) = 2
2
n=1 n
2
(20)
A continuacion, en la siguiente pagina, mostramos las soluciones numericas del Problema 1. para
diferentes valores de r 2 . El error lo definimos por la siguiente formula [4]
error =
(21)
La Diferencia se define por la resta entre el valor numerico de la funcion U (xi , yj ) obtenido por el
MDF y el valor analtico de la misma U (x, y).
Los resultados numericos de los problemas analizados y presentados en este trabajo son para los siguientes casos:
r = 0.48, r = 0.50, r = 0.30, r = 0.70. En las Tablas Dif significa Diferencia.
Caso 1.
1.r = 0.30, h = 0.1, k = 0.0030
x
t
u(xi , tj )
U(x, t)
Dif
0.100
0.006
0.2000
0.2000
0.0000
0.200
0.006
0.4000
0.3996
0.0004
0.300
0.006
0.6000
0.5941
0.0059
0.400
0.006
0.7640
0.7570
0.0070
0.500
0.006
0.8320
0.8252
0.0068
Error
0.0069
0.1024
0.9865
0.9218
0.8250
i
1
2
3
4
5
j
1
1
1
1
1
Caso 1.
2.r = 0.48, h = 0.1, k = 0.0048
x
t
u(xi , tj )
U(x, t)
Dif
0.100
0.005
0.2000
0.2000
0.0000
0.200
0.005
0.4000
0.3999
0.0001
0.300
0.005
0.6000
0.5970
0.0030
0.400
0.005
0.8000
0.7686
0.0314
0.500
0.005
0.8080
0.8436
0.0356
Error
0.0010
0.0300
0.4987
4.0841
4.2254
.
.
.
10
10
10
10
10
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.030
0.030
0.030
0.030
0.030
.
.
.
0.1829
0.3512
0.4890
0.5803
0.6123
.
.
.
0.1812
0.3484
0.4857
0.5770
0.6091
.
.
.
0.0016
0.0028
0.0033
0.0032
0.0032
.
.
.
0.9081
0.8046
0.6724
0.5621
0.5193
.
.
.
1
2
3
4
5
.
.
.
5
5
5
5
5
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.024
0.024
0.024
0.024
0.024
.
.
.
0.1898
0.3754
0.5093
0.6261
0.6389
.
.
.
0.1891
0.3658
0.5141
0.6146
0.6504
.
.
.
0.0007
0.0096
0.0048
0.0115
0.0115
.
.
.
0.3766
2.6170
0.9356
1.8742
1.7648
.
.
.
1
2
3
4
5
.
.
.
30
30
30
30
30
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.090
0.090
0.090
0.090
0.090
.
.
.
0.1033
0.1964
0.2704
0.3179
0.3342
.
.
.
0.1030
0.1960
0.2698
0.3171
0.3335
.
.
.
0.0002
0.0005
0.0006
0.0007
0.0008
.
.
.
0.2409
0.2376
0.2334
0.2301
0.2288
.
.
.
1
2
3
4
5
.
.
.
20
20
20
20
20
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.096
0.096
0.096
0.096
0.096
.
.
.
0.0969
0.1827
0.2536
0.2956
0.3135
.
.
.
0.0971
0.1847
0.2542
0.2989
0.3143
.
.
.
0.0002
0.0020
0.0006
0.0033
0.0008
.
.
.
0.2344
1.1016
0.2441
1.1113
0.2501
.
.
.
1
2
3
4
5
.
.
.
38
38
38
38
38
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.114
0.114
0.114
0.114
0.114
.
.
.
0.0814
0.1548
0.2130
0.2504
0.2633
.
.
.
0.0813
0.1547
0.2129
0.2502
0.2631
.
.
.
0.0001
0.0001
0.0002
0.0002
0.0002
.
.
.
0.0758
0.0752
0.0744
0.0738
0.0736
.
.
.
1
2
3
4
5
.
.
.
25
25
25
25
25
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.120
0.120
0.120
0.120
0.120
.
.
.
0.0756
0.1446
0.1980
0.2340
0.2447
.
.
.
0.0766
0.1458
0.2006
0.2359
0.2480
.
.
.
0.0010
0.0011
0.0027
0.0018
0.0033
.
.
.
1.3250
0.7712
1.3266
0.7728
1.3276
Error
0.0014
0.0382
0.5693
4.3469
4.8099
i
1
2
3
4
5
j
1
1
1
1
1
Caso
x
0.100
0.200
0.300
0.400
0.500
4.
j
1
1
1
1
1
Caso
x
0.100
0.200
0.300
0.400
0.500
3.
i
1
2
3
4
5
Error
0.022
0.212
1.491
7.023
11.24
.
.
.
1
2
3
4
5
.
.
.
6
6
6
6
6
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.030
0.030
0.030
0.030
0.030
.
.
.
0.1875
0.3438
0.5000
0.5625
0.6250
.
.
.
0.1812
0.3484
0.4857
0.5770
0.6091
.
.
.
0.0063
0.0046
0.0143
0.0145
0.0159
.
.
.
3.4570
1.3281
2.9338
2.5161
2.6058
.
.
.
1
2
3
4
5
.
.
.
5
5
5
5
5
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.035
0.035
0.035
0.035
0.035
.
.
.
0.1328
0.4576
0.2268
0.8898
0.1860
.
.
.
0.1741
0.3335
0.4630
0.5481
0.5778
.
.
.
0.0413
0.1242
0.2362
0.3417
0.3918
.
.
.
23.73
37.24
51.01
62.34
67.80
.
.
.
1
2
3
4
5
.
.
.
18
18
18
18
18
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.090
0.090
0.090
0.090
0.090
.
.
.
0.10490
0.18978
0.27466
0.30708
0.33951
.
.
.
0.1030
0.1960
0.2698
0.3171
0.3335
.
.
.
0.0019
0.0062
0.0049
0.0101
0.0060
.
.
.
1.8338
3.1551
1.8188
3.1715
1.8096
.
.
.
1
2
3
4
5
.
.
.
8
8
8
8
8
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.056
0.056
0.056
0.056
0.056
.
.
.
0.5529
-0.582
1.6757
-1.209
2.2456
.
.
.
0.1436
0.2735
0.3771
0.4439
0.4670
.
.
.
0.4093
0.8551
1.2985
1.6526
1.7786
.
.
.
285.00
312.58
344.31
372.26
380.82
.
.
.
1
2
3
4
5
.
.
.
24
24
24
24
24
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.120
0.120
0.120
0.120
0.120
.
.
.
0.07764
0.14044
0.20325
0.22724
0.25123
.
.
.
0.0766
0.1458
0.2006
0.2359
0.2480
.
.
.
0.0010
0.0053
0.0026
0.0086
0.0032
.
.
.
1.3101
3.6490
1.3085
3.6506
1.3075
.
.
.
1
2
3
4
5
.
.
.
11
11
11
11
11
.
.
.
0.100
0.200
0.300
0.400
0.500
.
.
.
0.077
0.077
0.077
0.077
0.077
.
.
.
-2.355
5.0544
-6.582
8.6683
-8.454
.
.
.
0.1171
0.2227
0.3067
0.3606
0.3793
.
.
.
2.4718
4.8317
6.8883
8.3077
8.8329
.
.
.
2111.4
2169.3
2246.2
2303.9
2329.4
i
1
2
3
4
5
j
2
2
2
2
2
.
.
.
1
2
3
4
5
t
0.005
0.005
0.005
0.005
0.005
t
0.007
0.007
0.007
0.007
0.007
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 1. r=0.30
2
Solucion aproximada
Solucion analitica
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 2. r=0.48
2
Solucion aproximada
Solucion analitica
1.5
U(x,t)
U(x,t)
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
x
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 3. r=0.50
0.8
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 4. r=0.70
2
Solucion aproximada
Solucion analitica
2
Solucion aproximada
Solucion analitica
1.5
1.5
U(x,t)
U(x,t)
0.6
x
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
Ahora analicemos un problema similar al Problema 1. pero empleando derivadas en las condiciones
de frontera
Problema 2.
U
2U
=
, 0 < x < 1, t > 0,
t
x2
U (x, 0) = 1, 0 x 1
= U, x = 0, t
U = U,
x
(22)
(23)
(24)
x = 1, t
Usemos el MDF en el problema y expresemos las condiciones de frontera mediante las diferencias
centrales, de tal manera que (22) queda as
ui,j+1 = ui,j + r(ui1,j 2ui,j + ui+1,j ).
(25)
(26)
(27)
(28)
Actuamos de la misma manera con la otra condicion del valor de frontera en x = 1, cuando i = 10
u10,j+1 = u10,j + r(u9,j 2u10,j + u11,j ),
(29)
u11,j u9,j
= u10,j ,
2h
de (30), tenemos
(30)
(31)
(32)
(33)
1
(34)
2
Hay que tener presente que tanto en el Problema 1. como en el Problema 2., con respecto a x = 12 ,
existe una simetra, la cual facilita, en este caso en particular, el calculo numerico.
tan =
Se ha demostrado que el esquema utilizado para el Problema 2. es convergente [3] siempre que
1
r 2+hx
.
En la siguiente pagina mostramos las Tablas que muestran los resultados numericos para el Problema
2.
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
1.
Error
0.11
0.70
0.02
0.00
0.00
0.00
i
0
1
2
3
4
5
j
1
1
1
1
1
1
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0002
0.0003
0.0008
0.0011
0.0013
0.0013
.
.
.
0.03
0.03
0.08
0.12
0.13
0.13
.
.
.
0
1
2
3
4
5
.
.
.
5
5
5
5
5
5
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0240
0.0240
0.0240
0.0240
0.0240
0.0240
.
.
.
0.8372
0.9227
0.9569
0.9887
0.9949
1.0000
.
.
.
0.8467
0.9172
0.9610
0.9841
0.9942
0.9968
.
.
.
0.0095
0.0055
0.0041
0.0046
0.0007
0.0032
.
.
.
1.12
0.60
0.43
0.46
0.07
0.32
.
.
.
0.7150
0.7800
0.8313
0.8683
0.8906
0.8981
.
.
.
0.0002
0.0000
0.0002
0.0004
0.0006
0.0006
.
.
.
0.03
0.00
0.03
0.05
0.06
0.07
.
.
.
0
1
2
3
4
5
.
.
.
21
21
21
21
21
21
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1008
0.1008
0.1008
0.1008
0.1008
0.1008
.
.
.
0.7144
0.7828
0.8317
0.8716
0.8917
0.9016
.
.
.
0.7165
0.7817
0.8331
0.8701
0.8924
0.8999
.
.
.
0.0021
0.0011
0.0014
0.0015
0.0007
0.0017
.
.
.
0.30
0.14
0.17
0.17
0.08
0.19
.
.
.
0.5999
0.6546
0.6981
0.7298
0.7490
0.7554
.
.
.
0.0004
0.0003
0.0002
0.0001
0.0001
0.0001
.
.
.
0.07
0.05
0.03
0.02
0.01
0.01
.
.
.
0
1
2
3
4
5
.
.
.
42
42
42
42
42
42
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2016
0.2016
0.2016
0.2016
0.2016
0.2016
.
.
.
0.6020
0.6563
0.7007
0.7320
0.7519
0.7578
.
.
.
0.6024
0.6573
0.7010
0.7328
0.7521
0.7585
.
.
.
0.0004
0.0010
0.0003
0.0008
0.0002
0.0007
.
.
.
0.06
0.15
0.04
0.11
0.03
0.10
i
0
1
2
3
4
5
j
1
1
1
1
1
1
.
.
.
0
1
2
3
4
5
.
.
.
10
10
10
10
10
10
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0300
0.0300
0.0300
0.0300
0.0300
0.0300
.
.
.
0.8308
0.9021
0.9496
0.9772
0.9906
0.9944
.
.
.
0.8311
0.9018
0.9488
0.9761
0.9893
0.9931
.
.
.
0
1
2
3
4
5
.
.
.
34
34
34
34
34
34
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1020
0.1020
0.1020
0.1020
0.1020
0.1020
.
.
.
0.7148
0.7800
0.8316
0.8688
0.8912
0.8987
.
.
.
0
1
2
3
4
5
.
.
.
68
68
68
68
68
68
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2040
0.2040
0.2040
0.2040
0.2040
0.2040
.
.
.
0.5995
0.6543
0.6979
0.7297
0.7489
0.7554
i
0
1
2
3
4
5
j
1
1
1
1
1
1
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0
1
2
3
4
5
.
.
.
4
4
4
4
4
4
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0200
0.0200
0.0200
0.0200
0.0200
0.0200
.
.
.
0.8691
0.9170
0.9775
0.9875
1.0000
1.0000
.
.
.
0.8585
0.9286
0.9695
0.9891
0.9967
0.9985
.
.
.
0
1
2
3
4
5
.
.
.
20
20
20
20
20
20
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1000
0.1000
0.1000
0.1000
0.1000
0.1000
.
.
.
0.7260
0.7741
0.8416
0.8646
0.9007
0.8952
.
.
.
0
1
2
3
4
5
.
.
.
41
41
41
41
41
41
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2050
0.2050
0.2050
0.2050
0.2050
0.2050
.
.
.
0.5868
0.6631
0.6868
0.7370
0.7385
0.7622
2.
Error
2.42
1.53
0.14
0.01
0.00
0.00
Error
2.70
1.62
0.16
0.01
0.00
0.00
i
0
1
2
3
4
5
j
1
1
1
1
1
1
Caso 4.
r = 0.70, h = 0.1, k = 0.0070
x
t
u(xi , tj )
U(x, t)
Dif
0.0000
0.0070
0.8600
0.9122
0.0522
0.1000
0.0070
1.0000
0.9752
0.0248
0.2000
0.0070
1.0000
0.9958
0.0042
0.3000
0.0070
1.0000
0.9996
0.0004
0.4000
0.0070
1.0000
1.0000
0.0000
0.5000
0.0070
1.0000
1.0000
0.0000
Error
5.72
2.55
0.43
0.04
0.00
0.00
.
.
.
0.0106
0.0116
0.0080
0.0016
0.0033
0.0015
.
.
.
1.24
1.24
0.83
0.16
0.33
0.15
.
.
.
0
1
2
3
4
5
.
.
.
3
3
3
3
3
3
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0210
0.0210
0.0210
0.0210
0.0210
0.0210
.
.
.
0.7576
0.9941
0.9314
1.0000
1.0000
1.0000
.
.
.
0.8554
0.9256
0.9673
0.9879
0.9961
0.9981
.
.
.
0.0978
0.0685
0.0359
0.0121
0.0039
0.0019
.
.
.
11.44
7.40
3.71
1.22
0.39
0.19
.
.
.
0.7176
0.7828
0.8342
0.8713
0.8936
0.9011
.
.
.
0.0085
0.0086
0.0074
0.0067
0.0071
0.0059
.
.
.
1.18
1.10
0.89
0.77
0.79
0.65
.
.
.
0
1
2
3
4
5
.
.
.
15
15
15
15
15
15
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1050
0.1050
0.1050
0.1050
0.1050
0.1050
.
.
.
-62.0875
56.4138
-46.6049
40.5358
-33.2094
32.9599
.
.
.
0.7113
0.7759
0.8270
0.8639
0.8862
0.8936
.
.
.
62.7987
55.6379
47.4319
39.6718
34.0956
32.0663
.
.
.
8829.33
7170.53
5735.09
4592.03
3847.51
3588.42
.
.
.
0.5989
0.6535
0.6970
0.7285
0.7477
0.7541
.
.
.
0.0121
0.0096
0.0102
0.0085
0.0092
0.0081
.
.
.
2.02
1.47
1.46
1.16
1.24
1.07
.
.
.
0
1
2
3
4
5
.
.
.
29
29
29
29
29
29
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2030
0.2030
0.2030
0.2030
0.2030
0.2030
.
.
.
-234318
212842
-195210
181872
-173483
170616
.
.
.
0.6009
0.6557
0.6993
0.7310
0.7503
0.7567
.
.
.
234319
212841
195211
181871
173484
170615
.
.
.
38992788
32458883
27913594
24878652
23123068
22546964
3.
10
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 1. r=0.30
2
Solucion aproximada
Solucion analitica
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 2. r=0.48
2
Solucion aproximada
Solucion analitica
1.5
U(x,t)
U(x,t)
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
x
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 3. r=0.50
0.8
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 4. r=0.70
2
Solucion aproximada
Solucion analitica
2
Solucion aproximada
Solucion analitica
1.5
1.5
U(x,t)
U(x,t)
0.6
x
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
11
El M
etodo Implicito Crank-Nicolson
Este metodo consiste en analizar la EDP en el punto (ih, (j + 21 )k) y aproximar la derivada con
respecto al tiempo por la diferencia central y la derivada con respecto a x por la media de las
diferencias centrales en los niveles de tiempo j, j + 1. Seg
un este metodo el problema (13) se expresa
en MDF as:
1 ui+1,j+1 2ui,j+1 + ui1,j+1 ui+1,j 2ui,j + ui1,j
ui,j+1 ui , j
.
(35)
=
+
k
2
h2
h2
Al operar (35) tenemos
rui1,j+1 + (2 + 2r)ui,j+1 rui+1,j+1 = rui1,j + (2 2r)ui,j + rui+1,j , r =
k
h2
(36)
Del lado izquierdo de (36) hay tres valores desconocidos y del lado derecho tres conocidos. Para
cada valor de t > 0, o j > 1, se forman N ecuaciones simultaneas con N incognitas, este sistema
puede ser resuelto , por ejemplo, por el Metodo de Gauss sin pivoteo.(Dibujar la celula)
Ahora resolvemos el Problema 1. por el metodo Crank-Nicolson; en la siguiente pagina mostramos
los resultados numericos para los distintos casos de r determinados en las anteriores tablas.
12
i
1
2
3
4
5
j
1
1
1
1
1
Caso
x
0.1000
0.2000
0.3000
0.4000
0.5000
Error
0.2354
0.4698
1.0978
1.5255
2.8773
i
1
2
3
4
5
j
1
1
1
1
1
Caso
x
0.1000
0.2000
0.3000
0.4000
0.5000
Error
0.3389
0.6499
1.2092
0.8688
1.6922
.
.
.
1
2
3
4
5
.
.
.
6
6
6
6
6
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0180
0.0180
0.0180
0.0180
0.0180
.
.
.
0.1770
0.3435
0.4874
0.5953
0.6549
.
.
.
0.1954
0.3816
0.5430
0.6561
0.6972
.
.
.
0.0183
0.0381
0.0557
0.0608
0.0423
.
.
.
9.3795
9.9848
10.2494
9.2634
6.0700
.
.
.
1
2
3
4
5
.
.
.
5
5
5
5
5
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0240
0.0240
0.0240
0.0240
0.0240
.
.
.
0.1756
0.3391
0.4770
0.5746
0.6194
.
.
.
0.1891
0.3658
0.5141
0.6146
0.6504
.
.
.
0.0135
0.0267
0.0371
0.0400
0.0310
.
.
.
7.1315
7.2965
7.2208
6.5085
4.7622
.
.
.
1
2
3
4
5
.
.
.
11
11
11
11
11
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0330
0.0330
0.0330
0.0330
0.0330
.
.
.
0.1410
0.2714
0.3809
0.4606
0.5037
.
.
.
0.1770
0.3394
0.4720
0.5594
0.5901
.
.
.
0.0359
0.0680
0.0911
0.0989
0.0864
.
.
.
20.3029
20.0451
19.3023
17.6717
14.6436
.
.
.
1
2
3
4
5
.
.
.
9
9
9
9
9
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0432
0.0432
0.0432
0.0432
0.0432
.
.
.
0.1389
0.2660
0.3700
0.4416
0.4739
.
.
.
0.1620
0.3092
0.4275
0.5044
0.5311
.
.
.
0.0231
0.0433
0.0575
0.0628
0.0572
.
.
.
14.2364
13.9879
13.4496
12.4537
10.7737
.
.
.
1
2
3
4
5
.
.
.
15
15
15
15
15
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0450
0.0450
0.0450
0.0450
0.0450
.
.
.
0.1156
0.2222
0.3113
0.3758
0.4106
.
.
.
0.1593
0.3040
0.4201
0.4954
0.5215
.
.
.
0.0437
0.0818
0.1088
0.1196
0.1109
.
.
.
27.4318
26.9190
25.9009
24.1373
21.2652
.
.
.
1
2
3
4
5
.
.
.
17
17
17
17
17
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0816
0.0816
0.0816
0.0816
0.0816
.
.
.
0.0832
0.1590
0.2207
0.2630
0.2820
.
.
.
0.1119
0.2129
0.2931
0.3446
0.3623
.
.
.
0.0287
0.0539
0.0723
0.0816
0.0804
.
.
.
25.6877
25.3273
24.6839
23.6810
22.1768
i
1
2
3
4
5
j
1
1
1
1
1
Caso
x
0.1000
0.2000
0.3000
0.4000
0.5000
Error
0.3488
0.6626
1.1920
0.7700
1.5667
i
1
2
3
4
5
j
1
1
1
1
1
Caso
x
0.1000
0.2000
0.3000
0.4000
0.5000
Error
0.4209
0.6749
0.7486
0.3877
0.4540
.
.
.
1
2
3
4
5
.
.
.
4
4
4
4
4
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0200
0.0200
0.0200
0.0200
0.0200
.
.
.
0.1836
0.3563
0.5042
0.6105
0.6590
.
.
.
0.1935
0.3766
0.5334
0.6418
0.6808
.
.
.
0.0099
0.0203
0.0291
0.0313
0.0218
.
.
.
5.1057
5.3891
5.4602
4.8731
3.2038
.
.
.
1
2
3
4
5
.
.
.
5
5
5
5
5
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0350
0.0350
0.0350
0.0350
0.0350
.
.
.
0.1662
0.3186
0.4433
0.5273
0.5604
.
.
.
0.1741
0.3335
0.4630
0.5481
0.5778
.
.
.
0.0079
0.0149
0.0196
0.0208
0.0175
.
.
.
4.5455
4.4577
4.2390
3.7969
3.0220
.
.
.
1
2
3
4
5
.
.
.
12
12
12
12
12
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0600
0.0600
0.0600
0.0600
0.0600
.
.
.
0.1130
0.2159
0.2997
0.3569
0.3822
.
.
.
0.1382
0.2631
0.3626
0.4267
0.4488
.
.
.
0.0252
0.0472
0.0628
0.0697
0.0665
.
.
.
18.2645
17.9370
17.3302
16.3473
14.8286
.
.
.
1
2
3
4
5
.
.
.
13
13
13
13
13
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0910
0.0910
0.0910
0.0910
0.0910
.
.
.
0.0895
0.1705
0.2356
0.2785
0.2952
.
.
.
0.1020
0.1940
0.2671
0.3140
0.3302
.
.
.
0.0125
0.0235
0.0315
0.0355
0.0350
.
.
.
12.2943
12.1169
11.8020
11.3149
10.5901
.
.
.
1
2
3
4
5
.
.
.
20
20
20
20
20
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1000
0.1000
0.1000
0.1000
0.1000
.
.
.
0.0666
0.1273
0.1767
0.2103
0.2252
.
.
.
0.0933
0.1776
0.2444
0.2873
0.3021
.
.
.
0.0267
0.0502
0.0677
0.0770
0.0769
.
.
.
28.6218
28.2939
27.7122
26.8120
25.4701
.
.
.
1
2
3
4
5
.
.
.
20
20
20
20
20
.
.
.
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1400
0.1400
0.1400
0.1400
0.1400
.
.
.
0.0513
0.0978
0.1352
0.1598
0.1694
.
.
.
0.0629
0.1197
0.1647
0.1936
0.2036
.
.
.
0.0116
0.0218
0.0295
0.0338
0.0342
.
.
.
18.3896
18.2221
17.9261
17.4705
16.7954
1.
3.
13
2.
4.
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 1. r=0.30
2
Solucion aproximada
Solucion analitica
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 2. r=0.48
2
Solucion aproximada
Solucion analitica
1.5
U(x,t)
U(x,t)
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
x
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 3. r=0.50
0.8
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 4. r=0.70
2
Solucion aproximada
Solucion analitica
2
Solucion aproximada
Solucion analitica
1.5
1.5
U(x,t)
U(x,t)
0.6
x
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
14
E Problema 2., puede ser resuelto aplicando el metodo Crank-Nicolson, teniendo en cuenta que
la expresion en MDF se ve as:
ui,j+1 ui,j
1 ui+1,j+1 2ui,j+1 + ui1,j+1 ui+1,j 2ui,j + ui1,j
(37)
=
+
k
2
h2
h2
Recordemos que expresamos la segunda derivada con respecto a x mediante la media de las derivadas
centrales, despues de algunas transformaciones tenemos:
rui1,j+1 + (2 + 2r)ui,j+1 rui+1,j+1 = rui1,j + (2 2r)ui,j + rui+1,j .
(38)
(39)
De (39) observamos que tenemos valores desconocidos como u1,j y u1,j+1 , para eliminarlos utilizamos las condiciones de frontera en x = 0 tanto para j como para j + 1 (tambien usando
diferencias centrales),
u1,j u1,j
= u0,j ,
(40)
2h
u1,j+1 u1,j+1
= u0,j .
(41)
2h
De lo anterior sigue
u1,j = u1,j 2hu0,j ,
(42)
y
u1,j+1 = u1,j+1 2hu0,j+1 .
(43)
Con (39) , (40) y (41) eliminamos los valores ficticios u1,j y u1,j+1 .
Por idea, algo parecido ocurre cuando x = 1, realizando un analisis similar podemos eliminar los
valores ficticios para i = N , en (39) tales como uN +1,j y uN +1,j+1 .
En la siguiente pagina mostramos los resultados numericos para el Problema 2. aplicando CrankNicolson.
15
Error
1.3376
0.1557
0.0399
0.0072
0.0009
0.0005
i
0
1
2
3
4
5
j
1
1
1
1
1
1
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
1.
j
1
1
1
1
1
1
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
1.
i
0
1
2
3
4
5
Error
0.7899
0.4042
0.0401
0.0249
0.0053
0.0027
.
.
.
0
1
2
3
4
5
.
.
.
5
5
5
5
5
5
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0150
0.0150
0.0150
0.0150
0.0150
0.0150
.
.
.
0.8794
0.9478
0.9814
0.9944
0.9984
0.9989
.
.
.
0.8755
0.9444
0.9802
0.9945
0.9988
0.9996
.
.
.
0.0039
0.0033
0.0012
0.0001
0.0004
0.0007
.
.
.
0.4480
0.3540
0.1220
0.0103
0.0424
0.0720
.
.
.
0
1
2
3
4
5
.
.
.
3
3
3
3
3
3
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0144
0.0144
0.0144
0.0144
0.0144
0.0144
.
.
.
0.8815
0.9498
0.9828
0.9950
0.9986
0.9991
.
.
.
0.8778
0.9465
0.9815
0.9951
0.9990
0.9997
.
.
.
0.0037
0.0033
0.0013
0.0000
0.0004
0.0006
.
.
.
0.4189
0.3504
0.1355
0.0046
0.0396
0.0557
.
.
.
0
1
2
3
4
5
.
.
.
34
34
34
34
34
34
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1020
0.1020
0.1020
0.1020
0.1020
0.1020
.
.
.
0.7130
0.7776
0.8278
0.8628
0.8822
0.8859
.
.
.
0.7150
0.7800
0.8313
0.8683
0.8906
0.8981
.
.
.
0.0021
0.0024
0.0035
0.0055
0.0084
0.0122
.
.
.
0.2887
0.3093
0.4258
0.6347
0.9409
1.3556
.
.
.
0
1
2
3
4
5
.
.
.
11
11
11
11
11
11
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0528
0.0528
0.0528
0.0528
0.0528
0.0528
.
.
.
0.7863
0.8562
0.9078
0.9421
0.9606
0.9650
.
.
.
0.7853
0.8551
0.9072
0.9425
0.9626
0.9691
.
.
.
0.0010
0.0011
0.0006
0.0004
0.0019
0.0040
.
.
.
0.1281
0.1295
0.0677
0.0443
0.2018
0.4148
.
.
.
0
1
2
3
4
5
.
.
.
67
67
67
67
67
67
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2010
0.2010
0.2010
0.2010
0.2010
0.2010
.
.
.
0.5948
0.6488
0.6911
0.7207
0.7373
0.7404
.
.
.
0.6030
0.6580
0.7017
0.7335
0.7528
0.7593
.
.
.
0.0082
0.0092
0.0107
0.0128
0.0156
0.0189
.
.
.
1.3648
1.3914
1.5215
1.7474
2.0694
2.4935
.
.
.
0
1
2
3
4
5
.
.
.
42
42
42
42
42
42
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2016
0.2016
0.2016
0.2016
0.2016
0.2016
.
.
.
0.5972
0.6515
0.6943
0.7248
0.7423
0.7466
.
.
.
0.6024
0.6573
0.7010
0.7328
0.7521
0.7585
.
.
.
0.0052
0.0057
0.0067
0.0080
0.0098
0.0119
.
.
.
0.8661
0.8746
0.9521
1.0933
1.2979
1.5698
Error
0.7294
0.4290
0.0369
0.0271
0.0061
0.0031
i
0
1
2
3
4
5
j
1
1
1
1
1
1
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
1.
j
1
1
1
1
1
1
Caso
x
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
1.
i
0
1
2
3
4
5
Error
0.1492
0.6379
0.0232
0.0457
0.0178
0.0098
.
.
.
0
1
2
3
4
5
.
.
.
3
3
3
3
3
3
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0150
0.0150
0.0150
0.0150
0.0150
0.0150
.
.
.
0.8791
0.9477
0.9816
0.9945
0.9984
0.9990
.
.
.
0.8755
0.9444
0.9802
0.9945
0.9988
0.9996
.
.
.
0.0035
0.0032
0.0014
0.0000
0.0004
0.0006
.
.
.
0.4030
0.3439
0.1421
0.0005
0.0409
0.0602
.
.
.
0
1
2
3
4
5
.
.
.
8
8
8
8
8
8
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.0560
0.0560
0.0560
0.0560
0.0560
0.0560
.
.
.
0.7808
0.8505
0.9026
0.9377
0.9574
0.9630
.
.
.
0.7799
0.8495
0.9018
0.9376
0.9582
0.9649
.
.
.
0.0009
0.0011
0.0008
0.0001
0.0008
0.0019
.
.
.
0.1164
0.1246
0.0851
0.0122
0.0829
0.1972
.
.
.
0
1
2
3
4
5
.
.
.
20
20
20
20
20
20
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1000
0.1000
0.1000
0.1000
0.1000
0.1000
.
.
.
0.7166
0.7817
0.8325
0.8684
0.8890
0.8941
.
.
.
0.7176
0.7828
0.8342
0.8713
0.8936
0.9011
.
.
.
0.0010
0.0011
0.0017
0.0029
0.0046
0.0069
.
.
.
0.1368
0.1405
0.2071
0.3323
0.5171
0.7674
.
.
.
0
1
2
3
4
5
.
.
.
15
15
15
15
15
15
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.1050
0.1050
0.1050
0.1050
0.1050
0.1050
.
.
.
0.7108
0.7755
0.8263
0.8626
0.8840
0.8902
.
.
.
0.7113
0.7759
0.8270
0.8639
0.8862
0.8936
.
.
.
0.0004
0.0004
0.0007
0.0013
0.0022
0.0034
.
.
.
0.0626
0.0542
0.0846
0.1492
0.2464
0.3776
.
.
.
0
1
2
3
4
5
.
.
.
40
40
40
40
40
40
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2000
0.2000
0.2000
0.2000
0.2000
0.2000
.
.
.
0.5991
0.6537
0.6967
0.7273
0.7450
0.7494
.
.
.
0.6040
0.6591
0.7029
0.7348
0.7541
0.7606
.
.
.
0.0049
0.0054
0.0062
0.0075
0.0092
0.0112
.
.
.
0.8077
0.8145
0.8870
1.0202
1.2137
1.4711
.
.
.
0
1
2
3
4
5
.
.
.
29
29
29
29
29
29
.
.
.
0.0000
0.1000
0.2000
0.3000
0.4000
0.5000
.
.
.
0.2030
0.2030
0.2030
0.2030
0.2030
0.2030
.
.
.
0.5984
0.6530
0.6962
0.7273
0.7457
0.7511
.
.
.
0.6009
0.6557
0.6993
0.7310
0.7503
0.7567
.
.
.
0.0025
0.0027
0.0031
0.0037
0.0045
0.0056
.
.
.
0.4217
0.4137
0.4436
0.5079
0.6059
0.7395
16
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 1. r=0.30
2
Solucion aproximada
Solucion analitica
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 2. r=0.48
2
Solucion aproximada
Solucion analitica
1.5
U(x,t)
U(x,t)
1.5
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
x
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 3. r=0.50
0.8
Cuadro de comparacion entre las soluciones aproximadas y analiticas por el MDF para el caso 4. r=0.70
2
Solucion aproximada
Solucion analitica
2
Solucion aproximada
Solucion analitica
1.5
1.5
U(x,t)
U(x,t)
0.6
x
0.5
0.5
0
0
0.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
17
Conclusi
on
En este artculo, Nosotros mostramos numericamente que para valores del parametro de estabilidad
r 12 , el metodo es estable(
util), en cambio para r > 12 el proceso es divergente. El metodo
implicito es convergente para cualquier r, teoricamente, sin embargo para r grandes, el metodo
genera soluciones oscilatorias.
Agradecimientos
Queremos agradecer al profesor Valdemir Garcia Ferreira por todos los consejos y sugerencias hechas
durante la realizacion de este trabajo. Extendemos nuestros agradecimientos tambien a la Universidad del Atlantico por el apoyo prestado.
Bibliografa
[1] N.N. Kalitkin, Chislennie Metodi, Nauka,Moskva,1978.
[2] G.U. Marchuk,Metodi vuichislitelnoi matematiki, Nauka, Moskva 1977.
[3] G.D. Smith, Numerical Solution of Partial DiffereNtial Equations.Finite Difference Methods.Third Edition.Oxford Apllied Mathematics and Computing Science Series.UK,2005.
` Solucao Numerica de Equacoes Diferenciais Parciais com
[4] V. Garcia Ferreira, Introducao A
Aplicacoes Em Dinamica dos Fluidos.Cursillo Universidad del Atlantico.2007.
Anexos
\\PROGRAMA QUE HALLA NUMERICAMENTE MEDIANTE MDF LA SOLUCION DEL PROBLEMA 1. POR EL METODO
EXPLICITO
#include <stdio>
#include<conio>
#include<iostream>
#include<math.h>
#define Pi 3.141592653589
using namespace std;
//r=0.48 explicito048.txt
//r=0.50 explicito050.txt
//r=0.70 explicito070.txt
//r=0.30 explicito030.txt
///funcion u(x,t) aproximada long double U(double x, double t)
{
long double z; int i; z=0;
for (i=1; i70;i++)
{
z = z+8/(Pi*Pi)*1/(i*i)*sin(Pi*i/2)*sin(i*Pi*x)*exp(-Pi*Pi*t*i*i);
}
18
return(z);
}
//fin funcion aproximada
////Main program
main ()
{
FILE *f;
f=fopen (explicito070.txt,w+);
int i,j,N,T,l,c;
float s1,s2,s3,r,h,k,x,xx,t,error,diferencia;
float zz;
float u[11][101];
r=0.70;
h=0.1;
k=r*h*h;
N=10;
T=50;
fprintf(f,r=%4.2f, h=%4.2f,k=%6.4f \n ,r,h,k);
fprintf(f,i j x
t
U[xi,yj] U[x,y] diferencia
error \n );
printf(i j x
t
U[xi,yj] U[x,y] diferencia
error \n );
for (l=0;l<=N+1;l++) //valores de U con las condiciones iniciales
{
x=l*h;
if ((x>=0 ) & (x<= 0.5))
{
u[l][0]=2*x;
fprintf(f,%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n ,l,0,x,0.0,u[l][0],u[l][0],0.0,0.0);
printf(%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,l,0,x,0.0,u[l][0],u[l][0],0.0,0.0);
}//fin del if
if ((x>0.5 ) & (x<=1))
{
u[l][0]=2*(1-x);
fprintf(f,%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n ,l,0,x,0.0,u[l][0],u[l][0],0.0,0.0);
printf(%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,l,0,x,0.0,u[l][0],u[l][0],0.0,0.0);
} //fin del if
} //fin del for
fprintf(f,%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,l,0,x,0.0,u[l][0],u[l][0],0.0,0.0);
printf(%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,l,0,x,0.0,u[l][0],u[l][0],0.0,0.0);
for (j=0;j=100;j++)
{
u[0][j]=u[10][j]=0;
t=(j+1)*k;
for (i=0;i<N+1;i++)
{
if ((i==0) || (i==10))
{
x=i*h;
fprintf(f,%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,i,j+1,x,t,u[i][j],u[i][j],0.0,0.0);
printf(%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,i,j+1,x,t,u[i][j],u[i][j],0.0,0.0);
}
else
{
x=i*h;
t=(j+1)*k;
s1=u[i-1][j];
s2=u[i][j];
s3=u[i+1][j];
u[i][j+1]=r*s1+(1-2*r)*s2+r*s3;
zz=U(x,t);
diferencia=fabs(zz-u[i][j+1]);
error=fabs((u[i][j+1]-zz)/zz)*100;
fprintf(f,%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,i,j+1,x,t,u[i][j+1],zz,diferencia,error);
printf(%d %d %5.3f %5.3f %6.4f %6.4f %6.4f %6.4f \n,i,j+1,x,t,u[i][j+1],zz,diferencia,error);
}
} // fin del for i
}/fin del for j
fclose(f);
printf(fin del programa );
19
20
\\PROGRAMA QUE HALLA NUMERICAMENTE MEDIANTE MDF LA SOLUCION DEL PROBLEMA 2. POR EL METODO
EXPLICITO
#include<stdio>
#include<conio>
#include<iostream>
#include<math.h>
#define Pi 3.141592653589
using namespace std;
// funcion original
double foriginal(double xx)
{
double t,zz;
t=cos(xx);
if (cos(xx)!=0)
{
zz= xx*tan(xx)-0.5;
return(zz);
}
else {
printf(Division por cero\n); exit(0); }
}
////
/// primera derivada de la solucion analitica
double fderivada(double xx )
{
double t,zz;
t= cos(xx);
if (t==0)
{
printf(hay un error\n); exit(1); }
else {
zz=tan(xx)+xx/(pow(cos(xx),2));
return(zz);
}
}
/// Metodo de Newton para hallar la soluci
on num
erica de xtanx=1/2
double Newton( int n )
{
int j;
double u1,u2,epsilon,zz,x1,x2;
double c,u3;
n=n-1;
c=(2*n+1)*Pi/2.0-Pi/180.0;
epsilon=0.001;
j=0;
x2=c;
do
{
x1=x2;
u2=foriginal(x1);
u3=fderivada(x1);
x2=x1-u2/u3;
j=j+1;
u1=fabs(x2-x1);
} while(u1epsilon);
zz=x2;
return(zz);
}//fin de la funcion newton///
//
long double U(double xx, double tt)
{
double alpha,z; int i;
z=0;
for (i=1;i<10;i++)
{
alpha=Newton(i);
z = 1/(cos(alpha)*(3+4*pow(alpha,2))) * exp(-4*pow(alpha,2)*tt)*cos(2*alpha*(xx-0.5))+z;
21
}
z=z*4;
return(z);
}
///////fin funcion analitica
///Main program
main ()
{
FILE *f;
f=fopen(explicitobvp030.txt,w+);
int i,j,N,T,l,m,c;
double s1,s2,s3,s4,s5,s6,s7,s8,r,h,k,x,xx,t,error,diferencia;
double zz,zzz;
float u[11][101];
r=0.30;
h=0.1;
k=r*h*h;
j=0;
N=10;
T=100;
fprintf(f,PARA EL CASO CUANDO BVP CON DERIVADAS PARCIALES\n);
fprintf(f,r=%4.2f, h=%6.4f,k=%6.k\n,r,h,k);
fprintf(f,i j x t U[i,j] Uanalitica diferencia error \n);
printf(i j x t U[i,j] Uanalitica diferencia error \n);
////vienen las condiciones iniciales
for (l=0; l<11;l++)
{
x=l*h;
s5=u[l][0]=1.0;
zzz=1.0;
diferencia=fabs(zzz-s5);
error=(diferencia)/zzz*100;
fprintf(f,%d %d %6.4f %6.4f %6.4f %6.4f %6.4f %4.2f \n, l,0,x,0.0,s5,zzz,diferencia,error);
printf(%d %d %6.4f %6.4f %6.4f %6.4f %6.4f %4.2f \n, l,0,x,0.0,s5,zzz,diferencia,error);
}
for (j=0;j<=T-1;j++)
{
///calculamos los valores en la frontera
s7=u[0][j];s8=u[1][j];
u[0][j+1]=u[0][j]+2*r*(u[1][j]-(1+h)*u[0][j]);
s5=u[0][j+1];
t=(j+1)*k;//ojo con la t
l=1;//un nuevo contador sirve cuando se vaya a buscar los valores simetricos
u[10][j+1]=u[10][j]+2*r*(u[9][j]-(1+h)*u[10][j]);
s6=u[10][j+1];
for (i=0;i<N+1;i++)
{
x=i*h;
if (i==0)
{
zzz=U(x,t);
diferencia=fabs(zzz-s5);
error=(diferencia)/zzz*100;
}
if (i==10)
{
zzz=U(x,t);
diferencia=fabs(zzz-s6);
error=(diferencia)/zzz*100;
}
if ((i<5) && (i!=0))
{
s1=u[i-1][j];
s2=u[i][j];
s3=u[i+1][j];
u[i][j+1]=s2+r*(s1-2*s2+s3);
s4=u[i][j+1];
}
22
if (i==5)
{
s1=u[i-1][j];
s2=u[i][j];
s3=u[i+1][j];
u[i][j+1]=s2+r*(s1-2*s2+s3);
s4=u[i][j+1];
}
if ((i>5) && (i!=N))
{
m=i-2*l;
u[i][j+1]=u[m][j+1];
l=l+1;
s4=u[i][j+1];
}
zzz=U(x,t);
diferencia=fabs(zzz-u[i][j+1]);
error=diferencia/zzz*100;
fprintf(f,%d %d %6.4f %6.4f %6.4f %6.4f %6.4f %4.2f \n, i,j+1,x,t,u[i][j+1],zzz,diferencia,error);
printf(%d %d %6.4f %6.4f %6.4f %6.4f %6.4f %4.2f \n, i,j+1,x,t,u[i][j+1],zzz,diferencia,error);
} // FIN DEL FOR PRINCIPAL DEL I
}//FIN DEL FOR PRINCIPAL DE J
fclose(f);
printf(fin del programa);
while((c = getchar()) != \n && c != EOF);
} ///fin del programa
23
\\PROGRAMA QUE HALLA NUMERICAMENTE MEDIANTE MDF LA SOLUCION DEL PROBLEMA 1. POR EL METODO
IMPLICITO
#include<iostream>
#include<stdlib>
#include<stdio>
#define Pi 3.141592653589
using namespace std;
///funcion u(x,t) aproximada
long double U(double x, double t)
{
long double z; int i;
z=0;
for (i=1; i<20;i++)
z = z+8/(Pi*Pi)*1/(i*i)*sin(Pi*i/2)*sin(i*Pi*x)*exp(-Pi*Pi*t*i*i);
return(z);
} ///////fin funcion aproximada //Main Program
int main()
{
int i,j,c,l,s;
FILE *f1;
int N=5;
float r=0.70;
float h=0.1;
float k=r*h*h;
float a[20][20],b[20][20],A[20],B[20],x,xx,diferencia,error,t,zzz;
f=fopen(implicito070.txt,w+);
fprintf(f1,PARA EL CASO IMPLICITO CRANCK NICHOLSON BVP \n);
fprintf(f1,r=%3.2f, h=%3.2f,k=%5.3f \n,r,h,k);
fprintf(f1, i j x t U[i,j] Uanalitica diferencia error \n);
printf(PARA EL CASO IMPLICITO CRANCK NICHOLSON BVP r=%3.2f, h=%3.2f\n,r,h);
printf( i j x t U[i,j] Uanalitica diferencia error \n);
for (i=0;i<=2*N+1;i++) {
x=i*h;
if ((x>=0 ) & (x<=0.5)) {
a[i][0]=2*x;
a[10-i][0]=a[i][0];
zzz=a[i][0];
diferencia=fabs(zzz-a[i][0]);
error=0;
fprintf(f1, %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,0,x,0.0,a[i][0],zzz,diferencia,error);
printf( %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,0,x,0.0,a[i][0],zzz,diferencia,error);
}
//fin del if
///////////////// por simetria tenemos
if ((x>0.5 ) & (x<1))
{
zzz=a[10-i][0];
diferencia=fabs(zzz-a[i][0]);
error=diferencia/zzz*100;
fprintf(f1, %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,0,x,0.0,a[i][0],zzz,diferencia,error);
printf( %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,0,x,0.0,a[i][0],zzz,diferencia,error);
}
//fin del if
if (x==1)
{
zzz=a[10-i][0];
diferencia=fabs(zzz-a[i][0]);
error=0;
fprintf(f1, %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,0,x,0.0,a[i][0],zzz,diferencia,error);
printf( %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,0,x,0.0,a[i][0],zzz,diferencia,error);
}
}
///fin del for i para llenarlos valores iniciales
for (j=0; j<20;j++) {
for (i=1;i<N+1;i++) {
24
if (i==1) {
b[i][i]=4.0;
b[i][i+1]=-1;
rellena de ceros el lado derecho
for (l=i+2;l<N+1;l++)
{
b[i][l]=0.0;
}
b[i][N+1]= (2-2*r)*a[i][j]+r*a[i+1][j];
printf(\n);
fprintf(f,\n);
}
if ((i>1) && (i!=N)) {
if(i>2) {
///rellena de ceros el lado izquierdo
for (l=1;l<i-1;l++)
{
b[i][l]=0.0;
}
}
b[i][i-1]=-1;
b[i][i]=4.0;
b[i][i+1]=-1;
for (l=i+2;l<N+1;l++)
{
b[i][l]=0.0;
}
b[i][N+1]= r*a[i-1][j]+(2-2*r)*a[i][j]+r*a[i+1][j];
}//fin del if i>1
if (i==N)
{
for (l=1;l<i-1;l++)
{
b[i][l]=0.0;
}
b[i][N-1]=-2.0*r;
b[i][N]=(2+2*r);
b[i][N+1]= r*a[i-1][j]+(2-2*r)*a[i][j]+r*a[i+1][j];
//a[i-1][j]+a[i+1][j];
printf(\n);
} }///fin del for I=1,N
printf(\n);
////METODO GAUSS SIN PIVOTEO-METOD PROGONKA//// //////// ///aplicamos progonka directa A[1]=-b[1][2]/b[1][1];
B[1]=b[1][N+1]/b[1][1];
for (i=2;i<N+1;i++)
{
if (i!= N)
{
A[i]=-b[i][i+1]/(b[i][i]+b[i][i-1]*A[i-1]);
}
B[i]=(b[i][N+1]-b[i][i-1]*B[i-1])/(b[i][i]+b[i][i-1]*A[i-1]);
}
//ahora hallamos las soluciones
//printf(\n);
xx=B[N];
//las soluciones son los elementos del vector a[i][j] a[N][j+1]=xx;
printf( a[%d][%d]= %6.4f \n,N,j+1,a[N][j+1]);
l=N-1;
while (l>0) {
s=10-l;
xx=A[l]*xx+B[l];
a[l][j+1]=xx;
a[s][j+1]=xx;
fprintf(f1, a[%d][%d]= %6.4f \n,l,j+1,a[s][j+1]);
printf( a[%d][%d]= %6.4f \n,l,j+1,a[l][j+1]);
l=l-1;
}
25
printf(\n);
t=k*(j+1);
for (i=0;i<2*N+1;i++) //ahora escribe las soluciones en el archivo de texto
{
x=i*h;
if ((i==0) || (i==2*N))
{
printf( %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,j+1,x,t,0.0,0.0,0.0,0.0);
fprintf(f1, %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,j+1,x,t,0.0,0.0,0.0,0.0);
}
else {
zzz=U(x,t);
diferencia=fabs(zzz-a[i][j+1]);
error=diferencia/zzz*100;
fprintf(f1, %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,j+1,x,t,a[i][j+1],zzz,diferencia,error);
printf( %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f \n, i,j+1,x,t,a[i][j+1],zzz,diferencia,error);
}
}
}//fin del j principal fclose(f);
fclose(f1);
while((c = getchar()) != \n && c != EOF);
}
\\PROGRAMA QUE HALLA NUMERICAMENTE MEDIANTE MDF LA SOLUCION DEL PROBLEMA 2. POR EL METODO
IMPLICITO
#include<iostream>
#include<stdlib>
#include<stdio>
#define Pi 3.141592653589
using namespace std;
//programa con el metodo Crank-Nicholson para bvp con derivadas double foriginal(double xx)
{
double t,zz;
t=cos(xx);
if (cos(xx)!=0)
{
zz= xx*tan(xx)-0.5;
return(zz);
}
else
{
printf(Division por cero \n);
exit(0);
}
}
////
/// primera derivada de la solucion analitica
double fderivada(double xx )
{
double t,zz;
t= cos(xx);
if (t==0)
{
printf(hay un error\n);
exit(1);
}
else
{
zz=tan(xx)+xx/(pow(cos(xx),2));
return(zz);
}
}
double Newton( int n )
{
int j;
double u1,u2,epsilon,zz,x1,x2;
double c,u3;
26
n=n-1;
c=(2*n+1)*Pi/2.0-Pi/180.0;
epsilon=0.001;
j=0;
x2=c;
do
{
x1=x2;
u2=foriginal(x1);
u3=fderivada(x1);
x2=x1-u2/u3;
j=j+1;
u1=fabs(x2-x1);
}
while(u1>epsilon);
zz=x2;
return(zz);
}
//fin de la funcion newton/// //la siguiente procedura halla el valor de la funcion analitica en un punto long double U(double xx, double
tt) {
double alpha,z; int i;
z=0;
for (i=1;i10;i++) {
alpha=Newton(i);
z = 1/(cos(alpha)*(3+4*pow(alpha,2))) * exp(-4*pow(alpha,2)*tt)*cos(2*alpha*(xx-0.5))+z;
}
z=z*4;
return(z);
}
///////fin funcion analitica
main ()
{
FILE *f;
f=fopen(imptobvpder070.txt,w+);
int i,j,N,T,l,m,c,s,ls,jj;
double s1,s2,s3,s4,s5,s6,r,h,k,x,xx,t,error,diferencia;
double zz,zzz;
float u[11][101],b[8][8],A[8],B[8];
r=0.70;
h=0.1;
k=r*h*h;
j=0;
N=5;
T=20;
fprintf(f,PARA EL CASO IMPLICITO CRANK-IMPLICITO CON BVP CON DERIVADAS PARCIALES\n);
fprintf(f,r=%5.3f, h=%5.3f,k=%6.4f\n,r,h,k);
fprintf(f, i j x t U[i,j] Uanalitica diferencia error \n);
printf( i j x t U[i,j] Uanalitica diferencia error \n);
////vienen las condiciones iniciales
for (l=0; l<11;l++)
{
x=l*h;
u[l][0]=1.0;
zzz=1.0;
diferencia=fabs(zzz-u[l][0]);
error=(diferencia)/zzz*100;
fprintf(f, %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %5.3f \n, l,0,x,0.0,u[l][0],zzz,diferencia,error);
printf( %d %d %6.4f %6.4f %6.4f %6.4f %6.4f %5.3f \n, l,0,x,0.0,u[l][0],zzz,diferencia,error);
}
//printf(\n);
for (j=0;j<100;j++)//ciclo del j principal
{
for (i=0;i<N+1;i++)
{
if (i==0)
{
b[i][i]=1+r+r*h;
b[i][i+1]=-r;
27
28
}
// printf(%6.4f, \n,u[10][50]);
/////////////fin progonka///////////////////// }
///fin del ciclo de j printf(fin del programa);
while((c = getchar()) != \n && c != EOF);
}
///fin del programa
29