How can I implement Newton raphson method in MATLAB for temperature modeling?

[ivanow]

how to write this equations in fortran or matlab

Energy(t) = Energy(t - dt) + (Solar_to_earth + IR_down - IR) * dt
INIT Energy = Initial_temperature*Heat_capacity
Solar_to_earth = Solar*(1-Atmospheric_absorption)
IR_down = IR_atmosphere
IR = Sigma*Temperature^4 * Seconds_per_year
Energy_in_atmosphere(t) = Energy_in_atmosphere(t - dt) + (IR + Solar_to_atmosphere - IR_down - IR_to_space) * dt
INIT Energy_in_atmosphere = Initial_temperature_atmosphere*Heat_capacity_atmosphere
IR = Sigma*Temperature^4 * Seconds_per_year
Solar_to_atmosphere = Solar*Atmospheric_absorption
IR_down = IR_atmosphere
IR_to_space = IR*(1-Emissivity) +IR_atmosphere
Albedo = IF(Albedo_feedback?=0) THEN 0.307168
ELSE Albedo_function
Albedo_feedback? = 1
Atmospheric_absorption = 0.1
Density_of_water = 1000 {kg/m^3}
Depth_of_mixed_layer = 50 {m}
E_0 = 0.75
Emissivity = if (H2O_feedback? = 0 ) then .856347
else
MAX(E_0 + .1 * T_1 * (1/T_0 - 1/Temperature),0)
H2O_feedback? = 1
Heat_capacity = Density_of_water*Depth_of_mixed_layer*Specific_heat_of_water
Heat_capacity_atmosphere = Mass_atmosphere*Specific_heat_atmosphere
Initial_temperature = 273.15 {K}
Initial_temperature_atmosphere = 273.15
IR_atmosphere = Sigma*Emissivity*T_atmosphere^4*Seconds_per_year
IR_down_norm = IR_down/Seconds_per_year
IR_norm = IR/Seconds_per_year
Mass_atmosphere = 101300/9.81 {kg/m^2}
Seconds_per_year = 365*24*60*60
Sigma = 5.6696e-8 {W/m^2/K^4}
Solar = Solar_constant/4*(1-Albedo)*Seconds_per_year
Solar_constant = IF(Stronger_sun?=0) THEN 1367 {W/m^2}
ELSE 1367 {W/m^2}*1.02
Solar_to_earth_norm = Solar_to_earth/Seconds_per_year
Specific_heat_atmosphere = 1004 {J/kg/K}
Specific_heat_of_water = 4218 {J/kg/K}
Stronger_sun? = 0
T_0 = 273.15 {K}
T_1 = 5423 {K}
T_atmosphere = Energy_in_atmosphere/Heat_capacity_atmosphere
T_Celsius = Temperature - 273.15
T_Celsius_atmosphere = T_atmosphere - 273.15
Temperature = Energy/Heat_capacity
Albedo_function = GRAPH(Temperature)
(200, 0.62), (201, 0.62), (202, 0.62), (203, 0.62), (204, 0.62), (205, 0.62), (206, 0.62), (207, 0.62), (208, 0.62), (209, 0.62), (210, 0.62), (211, 0.62), (212, 0.62), (213, 0.62), (214, 0.62), (215, 0.62), (216, 0.62), (217, 0.62), (218, 0.62), (219, 0.62), (220, 0.62), (221, 0.62), (222, 0.62), (223, 0.62), (224, 0.62), (225, 0.62), (226, 0.62), (227, 0.62), (228, 0.62), (229, 0.575), (230, 0.554), (231, 0.538), (232, 0.525), (233, 0.514), (234, 0.504), (235, 0.495), (236, 0.486), (237, 0.478), (238, 0.471), (239, 0.464), (240, 0.457), (241, 0.451), (242, 0.445), (243, 0.439), (244, 0.434), (245, 0.428), (246, 0.423), (247, 0.418), (248, 0.414), (249, 0.409), (250, 0.405), (251, 0.401), (252, 0.397), (253, 0.393), (254, 0.389), (255, 0.385), (256, 0.381), (257, 0.378), (258, 0.374), (259, 0.371), (260, 0.368), (261, 0.365), (262, 0.361), (263, 0.359), (264, 0.356), (265, 0.353), (266, 0.35), (267, 0.347), (268, 0.345), (269, 0.343), (270, 0.34), (271, 0.338), (272, 0.335), (273, 0.333), (274, 0.331), (275, 0.329), (276, 0.327), (277, 0.325), (278, 0.323), (279, 0.321), (280, 0.32), (281, 0.318), (282, 0.316), (283, 0.315), (284, 0.313), (285, 0.312), (286, 0.31), (287, 0.309), (288, 0.308), (289, 0.307), (290, 0.306), (291, 0.305), (292, 0.304), (293, 0.303), (294, 0.302), (295, 0.301), (296, 0.301), (297, 0.3), (298, 0.3), (299, 0.3), (300, 0.3), (301, 0.3), (302, 0.3), (303, 0.3), (304, 0.3), (305, 0.3), (306, 0.3), (307, 0.3), (308, 0.3), (309, 0.3), (310, 0.3), (311, 0.3), (312, 0.3), (313, 0.3), (314, 0.3), (315, 0.3), (316, 0.3), (317, 0.3), (318, 0.3), (319, 0.3), (320, 0.3), (321, 0.3), (322, 0.3), (323, 0.3), (324, 0.3), (325, 0.3), (326, 0.3), (327, 0.3), (328, 0.3), (329, 0.3), (330, 0.3), (331, 0.3), (332, 0.3), (333, 0.3), (334, 0.3), (335, 0.3), (336, 0.3), (337, 0.3), (338, 0.3), (339, 0.3), (340, 0.3), (341, 0.3), (342, 0.3), (343, 0.3), (344, 0.3), (345, 0.3), (346, 0.3), (347, 0.3), (348, 0.3), (349, 0.3), (350, 0.3)


time spec = 0 to 25 dt=0.01

the model calculate temperature for 11 values of initial temperatures

 

[gsal]

aren't YOU going to take a shot at it, first?
Or you expect one of us to do all the work for you?

'cause the second option is not quite the way things work around here.

So, please more specific as to what you are asking. Or better yet, first take a shot at it and THEN if you have problems there...we can help.

Also, please use the code - /code tags in square brackets to enclose source code to ease reading it.

 

[ivanow]

i write the entire code in MATLAB but i failed to calculed my temperature using different values of initial temperature. what tipe of iteration i need to use ?

 

[gsal]

Sorry, ivanow, I am afraid I don't know what to tell you...I don't even know what the original Stella code does...heck, I don't even know what Stella is.

But that's besides the point, I think you need to learn how to better ask a question in this forum.

You may have to start by posting your MATLAB code, indicating what the input is and what the expected output is and what it is that you are currently getting, etc...that could be a good start...starting to get the picture?

Other than that, I don't do matlab, I do fortran.

Good luck. Hopefully somebody else (matlab user) would come along.

 

[Bill Simpson]

Stella is the educational version of a system simulation tool.
http://www.iseesystems.com/softwares/Education/StellaSoftware.aspx

There are two related tasks to accomplish the task.

First you are sort of compiling one language into another, you need to know the rules about what are legal variable names, the syntax for assignments, for expressions, for if-then, and things like that. Most of the Stella code looks fairly straight forward to do this step.

Second you have to compile some statements at a higher level. Stella can numerically solve differential equations with feedback loops and inputs. You need to understand how to translate Stella code into Matlab's solving code. Most Stella diagrams don't really explicitly mention time, it is just implied and shows up on the graphs of the system behavior.

So, for starting any new project or language it is usually a good idea to start small. Can you write a VERY simple Stella diagram, possibly something where everything is constant and can you compile that into Matlab code that gives the same result? Then try adding an if-then and see if you can get that to translate correctly. Then try a very simple differential equation system and see if you can get that to work in Matlab. After you have gotten all that to work then you can throw all your equations into this and see if the result works and be able to figure out why not.

 

[ivanow]

i made it . i translated the entire code using differential equation solved whit a simple for where k was the size of time and if else iteration . now the problem is to implement Newton raphson method that in MATLAB is too difficult. can you give me some hints about Newton raphson method ? sorry for my english