mathematical modeling/differential equations

[miron1001]

Hello everyone.
I've just learned the existance of freebasic and now have to write a program of mathematical model of heat transfer in solar collector. But it doesnt matter a lot. What really matters is that i have all the equations, including differential ones. But dont know how to describe them in code.
Does everybody knows how to do it?
I would be happy if you give me some tips.
All the equations are below.

p*w*F=M=const
p=const
p*w*dw/dz=-dP-L*dz/l*p*w^2/2+p*q*sina*dz
p*w*F*d/dz*(i+w^2/2)=Qk
L*d^2*T/dz^2=Qp-Qk
M1*T1+M2*T(Z=l)=M*T(Z=0)

[Tourist Trap]

Hello everyone.
I've just learned the existance of freebasic and now have to write a program of mathematical model of heat transfer in solar collector. But it doesnt matter a lot. What really matters is that i have all the equations, including differential ones. But dont know how to describe them in code.
Does everybody knows how to do it?
I would be happy if you give me some tips.
All the equations are below.

p*w*F=M=const
p=const
p*w*dw/dz=-dP-L*dz/l*p*w^2/2+p*q*sina*dz
p*w*F*d/dz*(i+w^2/2)=Qk
L*d^2*T/dz^2=Qp-Qk
M1*T1+M2*T(Z=l)=M*T(Z=0)

Hello,

the problem is just a matter of time that I have not to discretize and put this into something that may run more or less. Sincerely, what you should do is find a program already made in basic or C, or mathematica, anything that would save us time, and we'll put the final touch to get it work in Freebasic. The equations given as you did are continuous, the step of discretization is wished because school memories are far :D

[Tourist Trap]

I would be happy if you give me some tips.

p*w*F=M=const
p=const
p*w*dw/dz=-dP-L*dz/l*p*w^2/2+p*q*sina*dz
p*w*F*d/dz*(i+w^2/2)=Qk
L*d^2*T/dz^2=Qp-Qk
M1*T1+M2*T(Z=l)=M*T(Z=0)

Ok for some tip anyway. Get rid of constant. They can be equalled to 1 in a first time.
So w= 1/F.
But second thing, which value controls the other? What is the function of time, or have you some other parameters varying independently like space coordinates? We dont see them here.

From what I see only w, P and z are implicitely designed to be parameters, and T only seems concerned by variations as far as M=cste

So maybe, make z only a variable to begin and do:

Code: Select all

function DeltaQ( DeltaZ as double) as double
     dim as double L, d, T, Qp, Qk
     Qp = 1
     return sqr(L*d^2*T/(Qp-Qk))
end function

dim as double dz
do
      dz += .1
      ? DeltaQ(dz)
loop until inkey()=chr(27)

Sorry but it's up to you now to bring something where one can see what are the function, of what, and how they should evolve (what range).

[lizard]

p*w*F=M=const
p=const
p*w*dw/dz=-dP-L*dz/l*p*w^2/2+p*q*sina*dz
p*w*F*d/dz*(i+w^2/2)=Qk
L*d^2*T/dz^2=Qp-Qk
M1*T1+M2*T(Z=l)=M*T(Z=0)

Basically you can write all kind of such equations directly in FreeBASIC code, fully described here:

https://www.freebasic.net/wiki/wikka.ph ... tPgOpIndex

[srvaldez]

I am no expert but it looks to me the equations would need to be solved first, then the symbolic solution (if any) could be evaluated with FB, if there's no closed-form solution then it would need to be evaluated using numerical methods.
I am simply too lazy to learn about differential equations in order to help.

[Tourist Trap]

I am no expert but it looks to me the equations would need to be solved first, then the symbolic solution (if any) could be evaluated with FB, if there's no closed-form solution then it would need to be evaluated using numerical methods.
I am simply too lazy to learn about differential equations in order to help.

This looks like a linear system (oups or maybe a degree 2 differential about q sighted) we can't say very much since we don't know the dependancies. Is it T(w), T(t), T(x,y,z)? And so on. Hard to solve anything if we don't have this information anyway.

[badidea]

I did not see this topic until 'Tourist Trap' responded. So this reply might be to late.
Anyway, it is difficult to help without more insight in the formula's and what each symbol represents e.g. "dw". Is "w" some width?
To show how you can do integration in (freebasic) code, I made the following:

Code: Select all

#include "crt/math.bi"
#include "string.bi"

'================ Simple demo of numerical integration of sine =================

dim as single angle, intSin 'use double for better accuracy
dim as single angleEnd = M_PI
dim as single angleStep = M_PI / 100

'integrate sine fron 0 to pi
for angle = 0 to angleEnd step angleStep
	intSin += sin(angle) * angleStep
next
print "Integral of sin(), 0 to pi:"; intSin

'See also: http://mathworld.wolfram.com/SineIntegral.html

'================ Lets integrate daily solar power on a panel ==================

'dayFraction = 0.0 to 1.0 = 00:00 to 24:00 hours (1 day)
function solarPowerDensity(dayFraction as single) as single
	dim as single maxPowerDensity = 1000 'W/m2 (surface perpendicalar)
	dim as single intensity = -cos(dayFraction * (M_PI *2))
	if intensity < 0 then intensity = 0 'sun light blocked by earth (night)
	return maxPowerDensity * intensity
end function

'A solar panel flat on earth surface
dim as single solarPanelSize = 1.0 * 1.5 'm2
dim as single solarPower, intSolarPower

for hour_ as integer = 0 to 24
	solarPower = solarPanelSize * solarPowerDensity(hour_ / 24)
	print str(hour_) + ":00 - Power [W]: " + format(solarPower, "0")
	intSolarPower += solarPower * 1 'multiple with 1 hour, use smaller step for better accuracy
next

print "Integrated solar power received [kWh]: " + format(intSolarPower / 1000, "0.0")

Assumptions: At the equator at September 21 or March 21, No clouds, etc.

[dodicat]

Examples of freebasic functions, numerically integrating (via simpson's rules') and a function root finder(bisection)

Code: Select all

screen 19
dim shared as integer xres,yres
screeninfo xres,yres

#macro sketch(_function,minx,maxx,miny,maxy,col)
For x As Double=minx To maxx Step (maxx-minx)/ordinates
    Dim As Double x1=(xres)*(x-minx)/(maxx-minx)
    Dim As Double y1=(yres)*(_function-maxy)/(miny-maxy)
   if x=minx then Pset(x1,y1),col else line -(x1,y1),col
Next x
#endmacro

#macro axis(colour)
Scope
    If Sgn(miny)<>Sgn(maxy) Then
        Line(0,(yres-(miny/(miny-maxy))*yres))-(xres,(yres-(miny/(miny-maxy))*yres)),colour 'x axis
        End if
        If Sgn(L)<>Sgn(U) Then
            Line(((L/(L-U))*xres),0)-(((L/(L-U))*xres),yres),colour 'y axis
            End if
                Draw String(0,yres/2),Str(L),colour
                Draw String(xres-8-8*(Len(Str(U))),yres/2),Str(U),colour
                Draw String(xres/2,0),Str(maxy),colour
                Draw String(xres/2,yres-16),Str(miny),colour
        End Scope
        #endmacro
        
'area by simpsons 
Function Area(fn As Function(x As Double) As Double,L As Double,U As Double,Ordinates As Integer=10000) As Double
    var n=Ordinates
    If n Mod 2=0 Then n=n+1 
    Var h=(U-L)/n
    Var Part1=fn(L+.5*h)
    Var Part2=0.0
    For i As Integer=1 To n-1
        Part1+=fn(L+h*i+h/2)
        Part2+=fn(L+h*i)
    Next i
    Function= (h/6)*(fn(L)+fn(U)+Part1*4+Part2*2) 
    'sketch part
    dim as double miny=1e100,maxy=-1e100
    For n As Double=L To U Step(U-L)/ordinates
        Dim As Double f=fn(n)
        If miny>f Then miny=f
        If maxy<f Then maxy=f
    Next n
    sketch(fn(x),L,U,miny,maxy,7)
    axis(6)
End Function

sub bisect( f as function(x as double) as double,min as double,max as double,byref O as double)
    dim as double last,st=(max-min)/100000,v
    for n as double=min to max step st
         v=f(n)
        if sgn(v)<>sgn(last) then print( n):O=n+st:exit sub
        last=v
    next
    end sub

sub roots( f as function(x as double) as double,min as double,max as double)
    dim as double last,O,st=(max-min)/10000000,v
    for n as double=min to max step st
         v=f(n)
        if sgn(v)<>sgn(last) and n>min then bisect(f,n-st,n,O):n=O
        last=v
    next
end sub

function norm(x as double) as double
    dim as double pi=4*atn(1),mu=0,sd=0,v=.05
    dim as double z=-(x-mu)^2/(2*v)
    return (1/sqr(2*pi*v)) * exp(z)
end function

function e(x as double) as double
    return exp(-x)
end function

function sine(x as double) as double
    return sin(x)
end function

function mixed(x as double) as double
    return norm(x)+e(x)+sine(x)
    end function



dim as double pi=4*atn(1)
locate 5
print "Area under a normal distribution =  ";Area(@norm,-1,1)
print "Press a key"
sleep
cls
locate 5
print "Area under sin(x) =";Area(@sine,0,pi)
print "press a key"
sleep
cls
locate 5
print "Area under exp(-x) = ";Area(@e,-2,10)
sleep
cls
locate 5
print "finding root of mixed"
roots(@mixed,-2,5)
'set to root of mixed
print "Area under mixed(x) up to an x intercept (root) = ";Area(@mixed,-2,3.1830625)
print "press a key to finish"
sleep